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EXECUTIVE  SUMMARY 


This  thesis  investigated  and  compared  alternative  signal  processing  techniques 
that  used  wavelet-based  methods  instead  of  traditional  frequency  domain  methods  for 
processing  measured  electromagnetic  pulse  (EMP)  waveforms.  The  work  was  based  on 
an  existing  procedure  developed  at  the  Naval  Air  Warfare  Center,  Aircraft  Division 
(NAWC-AD),  Patuxent  River,  MD  that  processed  signals  using  Fourier  transforms  and 
manual  preprocessing  techniques.  This  process  involved  the  capture,  preprocessing  and 
equalization  of  EMP  signals  that  penetrated  the  exterior  aircraft  shell. 

The  primary  focus  of  the  work  reported  in  this  thesis  was  equalization  and 
filtering  techniques  for  processing  EMP  signals  in  additive  white  noise.  The  physical 
model  of  the  NAWC-AD  test  procedure  consisted  of  receiver  and  digital  recording 
equipment  which  can  distort  the  desired  EMP  waveform’s  characteristics  prior  to  digital 
storage.  In  order  to  minimize  the  error  between  the  actual  and  recorded  signal,  channel 
equalization  and  noise  filtration  were  required.  This  signal  equalization  was  conducted  at 
the  sub-band  level  through  the  use  of  infinite  impulse  response  (HR)  filters  and  actual 
channel  response  characteristics.  Wavelet  based  decomposition  and  synthesis  provided 
the  mechanism  for  this  sub-band  analysis  to  occur.  Using  the  wavelet  techniques,  a  brief 
investigation  of  signal  de-noising  through  wavelet  thresholding  was  also  conducted.  This 
effort  led  to  the  development  of  a  number  of  different  equalization  techniques  for 
comparison  against  the  NAWC-AD  Fourier  equalization  technique. 

This  thesis  also  addressed  and  provided  viable  methods  for  signal  extraction  and 
DC  bias  removal  for  a  given  measured  EMP  waveform.  The  NAWC-AD  method 
employed  manual  procedures  to  identify  and  extract  the  recorded  wavefonn  based  on 
operator  visual  inspection.  The  use  of  time  averaged  energy  computation  with  an 
appropriately  selected  threshold  provided  the  basis  for  pre-processing  steps  with 
consistent  results.  By  automating  these  procedural  steps,  consistency  in  testing  results 
could  improve  the  overall  aircraft  testing  process.  The  mean  squared  error  is  used  as  the 
basis  for  the  comparison  of  the  effectiveness  of  the  equalization  algorithms. 
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It  was  found  that  wavelet  techniques  provided  results  that  were  as  good  or  better 
than  traditional  Fourier  techniques  employed  by  NAWC-AD.  In  systems  with  additive 
noise,  wavelet-based  techniques  exceeded  the  perfonnance  of  the  Fourier-based  methods 
and  surpassed  them  when  de-noising  techniques  were  used.  Results  suggest  that  the 
current  method  employed  by  NAWCAD  works  well  in  a  low  noise  (high  SNR) 
environment  based  on  the  channel  (equipment)  response  provided.  Significant 
improvement  in  performance  of  the  proposed  wavelet  methods  with  de-noising  is 
obtained  at  low  SNRs,  a  desirable  outcome  especially  when  the  aircraft  testing  is 
conducted  under  noisy  conditions. 


I.  INTRODUCTION 


The  Naval  Air  Warfare  Center,  Aircraft  Division  (NAWC-AD)  Patuxent  River, 
MD  conducts  testing  of  aircraft  against  the  effects  of  electromagnetic  pulse  (EMP) 
waveforms.  The  purpose  of  the  testing  is  to  detect  the  radio  frequencies  that  penetrate 
through  the  aircraft  shell  and  measure  EMP  waveforms  at  the  locations  of  critical  and 
sensitive  equipment.  It  is  then  determined  whether  this  measured  EMP  waveform  is 
within  safe  and  required  tolerances  for  the  aircraft.  In  efforts  to  upgrade  system 
components,  NAWC-AD  initiated  an  investigation  of  alternative  wavefonn  processing 
techniques  to  improve  the  effectiveness  of  this  type  of  aircraft  testing.  Present  waveform 
processing  techniques  involve  manual  processing  of  measured  EMP  wavefonns  to  extract 
the  signals,  equalization  of  signals  in  the  frequency  domain,  and  removal  of  the  DC  bias 
of  the  received  waveform. 

In  addition  to  manual  techniques  used  by  the  NAWC-AD,  the  signal  processing  is 
accomplished  through  Fourier  transforms,  which  exploit  the  frequency  characteristics  of 
a  signal.  Research  in  “Transform  Coded  eXcitation”  (TCX)  coding  has  shown  that  results 
for  “highly  non-stationary  signals,  such  as  percussions,”  are  not  optimal  using  discrete 
Fourier  transforms  [1].  Due  to  the  non-stationary  nature  of  the  EMP  waveforms,  it  is 
possible  that  alternative  transfonn  techniques,  specifically  wavelet  based-transforms,  may 
provide  better  waveform  processing  performance. 

A.  THESIS  OBJECTIVE 

This  thesis  presents  methods  for  improved  EMP  signal  processing  techniques  that 
would  eliminate  manual  processing  and  result  in  more  optimal  performance  through 
wavelet-based  transforms.  The  equalization  of  a  signal  through  the  use  of  discrete 
wavelet  transforms  provides  more  flexibility  through  the  ability  of  filtering  and 
processing  signals  at  the  sub-band  level  as  well  as  integration  of  accepted  de-noising 
algorithms. 

This  thesis  discusses  the  benefits  and  advantages  of  sub-band  processing, 

including  the  ability  to  allow  multiple  sampling  rates  and  de-noising.  Additionally,  it 

provides  a  possible  algorithm  for  acquisition  of  a  signal  using  time  average  energy 
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computations.  Performance  of  proposed  techniques  is  evaluated  through  MATLAB 
simulations  using  measured  EMP  waveforms  with  and  without  additive  noise. 


B.  THESIS  ORGANIZATION 

Chapter  II  discusses  the  signal  processing  background  required  to  address 
manipulation  and  processing  of  measured  EMP  waveforms.  Chapter  III  describes  the 
aircraft  testing  model  and  characteristics  of  the  given  measured  EMP  waveform  data  sets. 
Chapter  IV  outlines  the  currently  used  procedures  for  signal  processing  in  the  system  and 
introduces  proposed  wavelet-based  methods  to  process  the  EMP  waveforms.  Chapter  V 
presents  the  results  of  the  proposed  methods  and  their  comparison  to  the  Fourier  method. 
Chapter  VI  includes  conclusions  and  topics  for  further  research.  Appendix  A  contains  the 
MATLAB  code  used  to  process  the  EMP  waveforms. 
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II.  SIGNAL  PROCESSING:  TECHNIQUES  AND  DOMAINS 


This  chapter  provides  the  background  for  signal  representations  in  the  time 
domain  and  frequency  domain.  It  provides  an  overview  of  the  Fourier  transform 
techniques  that  are  the  foundation  of  the  currently  used  algorithms.  It  also  introduces  the 
wavelet  transforms  and  basic  de-noising  techniques.  Finally,  this  chapter  discusses  the 
basic  concepts  of  channel  filtering  and  equalization. 

A.  SIGNAL  DOMAINS 

A  given  signal,  in  particular  an  electro-magnetic  waveform,  can  be  analyzed  by 
using  a  number  of  different  methods.  The  use  of  mathematical  transformations  can  lead 
to  more  versatile  representation  of  signals  for  processing.  For  example,  although  some 
characteristics  of  signals  may  not  be  evident  in  the  time  domain,  definitive  signal 
properties  may  be  clearly  observed  in  the  frequency  domain.  Fourier  transform  and 
wavelet  decomposition  are  two  methods  of  expressing  time  domain  signals  in  another 
domain  in  order  to  extract  information  as  well  as  perform  processing  with  greater  ease  in 
a  computationally  efficient  manner. 


1.  Fourier  Transform 

Using  the  Fourier  transform,  it  is  possible  to  represent  a  given  time-domain  signal 
as  a  function  of  all  frequencies  present  in  it.  This  is  accomplished  by  representing  a 
signal  x(t)  as  a  linear  combination  of  complex  exponentials  [2],  given  by 


X{co)=  |  x(t)e-jco‘dt  (2.1) 

-oo 

where  X(oj)  is  the  transformed  signal  and  co  is  the  radian  frequency. 

An  inverse  transform  can  be  applied  to  transform  a  frequency  domain  signal  X(co) 
to  a  time-domain  signal  x(t),  given  by 


a  +co 

x(t)  =  —  f  X(0J)ejmtd0J. 


(2.2) 
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The  continuous-time  Fourier  transform,  however,  can  only  be  applied  to  signals 
of  infinite  duration  and  continuous  time.  Due  to  the  wide-spread  use  of  digital 
technologies  in  modern  signal  processing  applications,  signals  are  commonly  acquired 
and  stored  digitally  as  a  set  of  data  points  or  a  vector.  For  this  reason,  it  is  necessary  to 
utilize  the  discrete  Fourier  transfonn. 


2.  Discrete  Fourier  Transform  (DFT) 

The  DFT  is  an  extension  of  the  Fourier  transfonn  to  process  digitized  signal 
sequences  of  finite  duration  and  provides  a  one-to-one  correspondence  of  a  time  domain 
signal  to  its  frequency  domain  counterpart  for  a  given  sampling  rate.  The  transform  also 
represents  the  time-domain  signal  as  a  linear  combination  of  complex  exponentials,  but 
these  exponentials  occur  at  discrete  frequency  locations  as  opposed  to  over  the 
continuous  spectrum. 

In  discrete  signal  processing,  input  signals  are  typically  assembled  as  data  vectors 
of  length  N.  If  Ts  is  the  sampling  period  of  the  digitized  signals,  the  duration  of  a  vector 
containing  N  samples  is  NTS.  The  DFT  of  N  samples  of  a  discrete-time  signal,  x\n\,  is 
given  by  [2] 

1  N~  1 

1 M  =  77  Z  ■ *  M  e-M2*/N)n ,  k  =  0, 1, . . . ,  N  - 1  (2.3) 

N  n=0 


where  n  and  k  represent  the  time  and  frequency  indices.  The  DFT  represents  the  signal  in 

2 Tck 

the  frequency  domain  as  a  function  of  the  digital  frequency,  cod  =  — ;  0  <  cod  <  In .  The 
digital  frequency  can  be  expressed  as  cod  =  where /is  the  analog  frequency  and 


1  F 

Fs  =  —  is  the  sampling  frequency.  The  width  of  each  DFT  frequency  bin  is  Af  =  . 


The  fast  Fourier  transform  (FFT)  is  a  computationally  efficient  algorithm  to 
implement  the  DFT.  The  FFT  utilizes  a  “divide  and  conquer”  technique  to  perfonn  the 
DFT  calculations  quickly  [3]  by  breaking  a  signal  sequence  down  systematically  into 
smaller  sequences.  For  the  FFT  to  be  effective,  the  sequence  length  must  be  an  integer 
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power  of  2.  This  can  quickly  be  achieved  by  appending  zeros  to  the  end  of  a  signal  vector 
called  “zero  padding,”  which  has  no  effect  on  the  characteristics  of  the  signal  being 
analyzed. 

Figure  1  shows  a  chirp  signal  in  the  time-domain  and  its  frequency-domain 
representation  (magnitude  only)  obtained  by  using  a  1 ,024-point  FFT  algorithm.  From 
Figure  1,  it  can  be  seen  that  the  signal  contains  stronger  low  frequency  components  than 
high  frequency  components.  However,  where  in  time  these  components  occur  cannot  be 
determined  from  the  frequency  response.  Other  techniques,  such  as  wavelet  transforms, 
not  only  provide  frequency  information  of  a  signal,  but  also  provide  corresponding  time 
information  as  well  when  transformed.  For  this  reason,  the  wavelet  transforms  are 
described  in  the  following  sections. 
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Figure  1.  FFT  of  Signal  with  additive  noise,  (a)  Original  Signal:  1,024-Point  Frequency 
Chirp  with  Additive  White  Noise,  (b)  FFT  (magnitude  response)  of  the  Chirp 

Signal. 
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3. 


Continuous-Time  Wavelet  Transform  (CWT) 


The  wavelet  transform  is  another  method  for  representing  signals.  The  wavelet 
transform  has  the  ability  to  represent  the  time-frequency  variations  of  signal  components 
in  multiple  discrete  frequency  bands,  which  can  improve  analysis  of  a  given  signal.  The 
CWT  of  a  signal  x{t)  is  defined  as  [4] 


+co 

W (w,s)  =  J  x(t)— j 


V 


t-u 


dt 


(2.4) 


where  y/(t)  is  the  wavelet  basis  function;  the  variable  u  represents  the  time  shifting  of  the 
basis  function  and  the  variable  s  represents  the  scaling  value  of  the  wavelet  basis 
function.  The  wavelets  basis,  y/(t),  has  a  time  average  of  zero: 

+oo 

=  0  (2.5) 


For  each  basis  function,  a  scaling  function,  (j){t ) ,  exists  such  that 


-rw 

|^(0|2  =  J  ItK-sOl 


c is 


(2.6) 


and  s  represents  the  scaling  value  of  the  function  for  radian  frequency  co  as  defined 
above.  The  function,  i//(t)  and  </>(t) ,  are  orthogonal  such  that 


J 


(t)dt  =  0 


(2.7) 


ensuring  that  their  projections  into  the  wavelet  domain  are  unique  and  consistent 
transformation  of  the  signal  without  loss  or  redundancy  of  information  [4]. 

This  transformation  is  more  versatile  than  the  Fourier  transfonn.  By  manipulating 
the  shifting  and  scaling  parameters,  a  signal  transfonned  by  wavelets  provides  a  multi¬ 
band  representation  that  can  detect  and  characterize  transient  components  with  a  zooming 
procedure  across  scales  [4].  Wavelet  types  and  their  advantages  are  discussed  further  in 
Chapter  IV. 
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As  with  the  Fourier  transform,  the  discrete  representation  of  a  signal  utilizing  the 
wavelet  transform  can  also  be  found  using  the  discrete  wavelet  transform.  Additionally, 
there  is  further  discussion  on  the  benefits  of  using  wavelet  transforms  over  Fourier 
transforms  in  the  following  sections. 

4.  Discrete  Wavelet  Transform  (DWT) 

The  process  of  wavelet  analysis  is  a  multi-step  procedure  in  which  a  signal  is 
decomposed  into  several  sub-band  components.  This  transform  can  be  written  as  a 
convolution  for  a  discrete  signal  x[n\  with  a  wavelet  y/ j[n\  as  follows: 

JV-l 

W^n,aJJ  =  y^x\m}//  ■  \m-n\ .  (2.8) 

m= 0 

In  implementation,  the  discrete  wavelet  transform  can  be  realized  using  a  series  of 
filter  pairs  that  decompose  a  signal  into  a  number  of  sub-band  signals  as  required.  A 
discrete-time  signal  is  applied  to  a  low-pass  filter  having  an  impulse  response  h'  and  a 
high-pass  filter  with  an  impulse  response  g' .  These  filters  complement  each  other  and 
divide  the  spectral  range  of  the  input  signal  x[«]  into  a  low  frequency  band  and  a  high- 
frequency  band,  respectively.  The  impulse  responses  of  these  filters  must  satisfy  certain 
properties  [4].  Let  h'[n\  and  g'[n ]  be  the  reversals  of  impulse  responses  h[n]  and  g[n\, 
respectively,  i.e., 

h\n\  =  h[-n] 

gVA  =  g[-»] 

which  in  turn  are  related  to  each  other  as  given  by  g[/?]  =  (—  l)1  "  A[1  —  w] .  The  reader  is 
referred  to  [4]  for  an  in-depth  discussion  of  the  properties. 

The  wavelet  and  scaling  function  of  the  DWT  can  be  recursively  generated  using 
these  filters  as  follows: 

+oo 

0j(n)  =  V2  X  h[pWj+i(2n~P) 

P=^°  (2.9) 

+oo  v  7 

¥}  («)  =  V2  X  g\.P¥J+ 1  (2»  -  P) 
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where  n  represents  the  time  indices  of  the  wavelet  function,  j  represents  the  resolution 
level,  and  the  term  2  in  <j)(2n  -  p)  indicates  a  decimation  operation  by  2. 


Let  ej[n\  and  fj[n\  represent  the  approximation  and  detail  coefficients, 
respectively,  which  are  obtained  as  the  following  inner  products  [4] 


ej[n]  =  (x[n\,</>j[n]) 
fj[n]  =  (x[n\,i//j[n]) 


(2.10) 


These  inner  products  in  turn  yield  recursive  operations  for  computing  the  approximation 
and  detail  coefficients  [4]: 


ej[n]  =  X  hlP~2n]ej-ilP~\ 

p=— oo 

+CO 

fj[n]=  X  g[P~2n^ej-i[P~\ 


where  j  defines  the  level  of  decomposition.  Figure  2  shows  a  schematic  of  the 
implementation  of  Equation  (2.1 1).  The  decimation  is  represented  by  ^  2  in  the  figure. 
Each  level  of  decomposition  results  in  approximation  and  detail  coefficients  denoted  ej 
and  fh  respectively.  Through  repetitive  iterations,  this  process  can  continue  starting  from 
the  original  signal,  e0,  down  to  where  f,  and  en  are  a  length  of  one. 


Figure  2.  Two-level  Wavelet  Decomposition  using  a  Tree-Structured  Filter  Rank  (After 

Ref.  [4]). 
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Conversely,  wavelet  reconstruction  is  the  process  of  taking  sub-band  coefficients 
and  rebuilding  the  signal  to  its  original  form.  The  signal  reconstruction  is  given  by  [4] 

+oo  +00 

ej[n\  =  X  p]ej+l[p]  +  Yj  g\-n~2PUj+i[p]  •  (2-12) 

P=— 00  P=— 00 


As  with  wavelet  decomposition,  reconstruction  employs  a  similar  tree-like  filter 
bank  structure  as  shown  in  Figure  3.  In  this  diagram,  the  filters  used  are  h  and  g,  the 
reversals  of  h '  and  g',  respectively.  The  coefficients  are  filtered  and  then  up-sampled  by 
two,  denoted  by  the  symbol  T  2 . 


fj+M 


e'J+2[n] 


e’,[n] 


Figure  3.  Wavelet  Reconstruction  (After  Ref.  [4]). 


By  using  different  wavelet  basis  functions  or  the  number  of  levels,  there  are 
multiple  possible  representations  of  time-domain  signal.  Figure  4  illustrates  the  results  of 
wavelet  decomposition  of  an  arbitrary  1 ,024-point  frequency  chirp  signal  with  additive 
noise  and  three  levels  of  decomposition  using  the  Haar-wavelet.  As  can  be  seen,  the  <?s 
coefficients  shown  in  Figure  4  (e)  contain  the  low  frequency  portion  of  the  signal  and  are 
of  length  128  points.  The  f\  coefficients  shown  in  Figure  4  (b)  are  the  high  frequency 
coefficients  of  length  512. 
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Magnitude  Magnitude  Magnitude 


Figure  4.  Detail  and  Approximation  Coefficients  of  a  Sinusoidal  Chirp  in  Additive  Noise 
Using  and  Haar- wavelet  Transform:  (a)  Given  Chirp  Signal  of  1,024  Points,  (b) 
Detail  Level  1  Coefficient,  f\.  (c)  Detail  Level  2  Coefficients,^-  (d)  Detail  level  3 
Coefficients,^,  (e)  Approximation  Level  3  Coefficients  e$. 


Reconstruction  can  occur  up  to  the  original  signal  or  can  occur  for  specific 
branches.  Selective  reconstruction  is  when  approximation  or  detail  coefficients,  e  j[n]  or 

fj\n] ,  are  reconstructed  as  in  Figure  3  while  all  other  sub-band  coefficients  are  set  to 

zero.  The  signal  is  reconstructed  through  h  and#  until  it  is  e"o[n].  The  signal  e"o[n]  is  the 
time-domain  representation  of  the  desired  sub-band. 

In  summary,  wavelet  decomposition  divides  the  frequency  domain  and  divides  it 
into  equal  sub-bands  successively,  i.e.,  it  decomposes  the  signal  into  a  high  frequency 
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component  and  a  low  frequency  component.  As  evident  in  Figure  4,  each  time  the  signal 
decomposes  to  the  next  level,  there  are  half  as  many  points  to  represent  the  data  than  at 
the  level  before.  The  wavelet  decomposition  in  Figure  4  gives  an  alternative 
representation  of  the  sinusoidal  chirp  in  additive  noise  compared  to  that  given  by  the  FFT 
in  Figure  1  (b).  The  fact  that  a  different  representation  of  the  same  signal  exists  leads  to 
the  possibility  of  exploiting  these  representations  using  techniques  not  available  to  the 
FFT-based  signal  processing. 


5.  Types  of  Wavelets 

Wavelets  have  greater  versatility  than  Fourier  transforms  because  they  allow  for 
multiple  types  of  basis  functions  and  multiple  levels  of  decomposition.  As  discussed, 
wavelet  basis  functions  are  associated  in  pairs  as  defined  in  Equations  (2.4)  and  (2.5), 
y/(t),  and  a  scaling  function  (f>  (/),  and  satisfy  the  relationship  in  Equation  (2.6).  As  long 
as  the  wavelet  has  a  time  average  of  zero  as  defined  in  Equation  (2.5)  and  the  scaling 
function  is  orthogonal  to  the  basis  wavelet  and  satisfies  Equation  (2.8)  and  (2.9),  it 
constitutes  a  valid  wavelet  pair  [5].  This  allows  virtually  no  limit  on  the  number  of 
wavelet  bases  available.  However,  in  the  literature,  some  popular  wavelet  families  have 
emerged  due  to  their  superior  performance. 

For  a  given  signal,  some  wavelets  produce  better  results  than  other  wavelets.  For 
example,  the  “Haar”  wavelet  is  a  square  wave  function,  as  depicted  in  Figure  5  (a)  and 
defined  as 


V(0  = 


fl  0  <  t  <  0.5 
[-1  0.5  <  f  <  1 


(2.13) 


with  the  ability  to  represent  a  square  wave  easily  and  accurately  using  only  a  few  signal 
samples.  The  corresponding  Haar  scaling  function  is  given  by 

fl  0  <  r  <  1 

(/>(t)  = 

[0  otherwise 


However,  for  a  typical  sinusoidal  signal,  the  “Haar”  may  not  be  the  best  basis  function  to 
utilize. 
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Wavelets  are  generally  classified  into  families.  Some  typical  families  are  Haar, 
Daubechie’s,  Coiflet,  and  Meyer.  With  exception  of  the  Haar  and  Meyer  wavelets  of  this 
sub  set,  the  other  wavelets  do  not  have  explicit  expressions  that  define  them.  Figure  5 
shows  a  plot  of  each  respective  wavelet  type.  Unless  sufficient  a  priori  information  is 
available  for  the  signal  and  its  correlation  to  a  wavelet,  it  cannot  be  determined  which 
wavelet  basis  function  will  have  the  best  results.  Typically,  wavelets  are  optimized  when 
there  is  a  minimum  number  of  wavelet  coefficients  that  have  values  that  are  much  greater 
than  zero  [4,  5]. 
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Figure  5.  Types  of  Wavelet  Basis  Functions:  (a)  Haar  Wavelet,  (b)  Daubechie’s  Order  3 

Wavelet,  (c)  Meyer  Wavelet,  (d)  Coiflet  Wavelet. 
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In  this  work,  a  subset  of  the  different  types  of  wavelet  basis  function  are  used  and 
compared.  Those  used  in  this  work  include  Haar,  Daubechie’s  order  3  and  15,  Coiflet, 
and  Meyer.  No  research  was  conducted  to  determine  the  best  wavelet  for  this  application. 


6.  De-Noising 


The  removal  of  noise  in  signals  can  be  accomplished  by  a  number  of  different 
ways  using  digital  or  analog  filters  that  eliminate  frequencies  outside  the  range  of  interest 
of  a  given  signal.  Effective  filters  can  be  designed  for  processing  signals  affected  by 
white  noise  with  sufficient  a  priori  knowledge  of  the  desired  signal  and  added  cost  of 
complexity.  For  wavelets,  a  simple  solution  that  does  not  require  knowledge  of  the  signal 
for  the  removal  of  noise  can  be  implemented  using  de -noising  or  thresholding  techniques. 
For  a  given  normalized  signal  vector,  x\ri\,  a  sample-by-sample  manipulation  called 
thresholding  is  applied.  Hard  thresholding  is  implemented  as 


UM  (f\4n]\  >  T 
{o  if\x[n]\<T 


and  for  soft  thresholding 


xST[n]  =  \ 


x\n\  -  T 
x\n\+T 


1° 


if  x[n]  >  T 
if  x[n\  <  T 
if\x[n]\<T 


(2.14) 


(2.15) 


where  T  is  an  arbitrarily,  empirically,  or  theoretically  detennined  threshold  value. 
Combining  Stein’s  Unbiased  Estimate  of  Risk  (SURE)  method  and  fixed  thresholding 

<jyJ 2  log  N  where  N  is  the  length  of  the  signal  vector  and  o  is  a  scaling  multiple,  is  a 

MATLAB  function  called  heuristic  de-noising  which  takes  advantage  of  both  methods. 
The  results  in  Chapter  V  used  this  heuristic  method;  the  SURE  method  is  used  unless  a 
high  signal  to  noise  ratio  is  detected  by  the  function  [6]. 

Figure  6  depicts  a  graphical  representation  of  the  transfer  characteristics  for  a 
signal’s  magnitude  over  a  normalized  scale. 
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Normalized  Threshold  value 


Normalized  Signal  x[n] 
Magnitude  Value 

(a) 


Normalized  Signal  x[«] 
Magnitude  Value 

(b) 


Normalized  Signal  x[n] 
Magnitude  Value 

(c) 


Figure  6.  Normalized  Signal  Values  versus  Hard  and  Soft  Thresholding  with  a  Threshold 
Value  of  Approximately  0.5  (After  Ref.  [5]).  (a)  Original  Signal  Values,  x[n] 
with  no  thresholding  or  T—  0.  (b)  Hard  Thresholded  Signal  Value,  xm[n\  for  T= 
0.5.  (c)  Soft  Thresholded  Signal  Value,  xst[ii\  for  T=  0.5. 


To  implement  de-noising,  a  signal  is  decomposed  into  approximation  and  detail 
coefficients  using  the  wavelet  transfonn.  These  resultant  coefficients  are  then  nonnalized 

and  a  threshold  value  is  calculated  using  the  SURE  method  or  ^2  log  N  to  minimize  the 

estimation  error  between  the  actual  noise  signal  and  the  signal  estimate  without  noise  [4]. 
Wavelet  coefficients  are  then  compared  against  this  threshold  and  the  applicable 
threshold  is  applied.  All  values  above  the  threshold  are  left  alone  in  the  case  of  hard 
thresholding  shown  in  Figure  6  (b);  in  the  case  of  soft  thresholding,  the  values  are  scaled 
lower  as  shown  in  Figure  6  (c).  The  assumption  behind  this  method  of  processing  is  that 
unless  the  magnitude  of  the  coefficient  exceeds  a  threshold  value,  it  is  a  contribution 
from  noise  and  should  be  removed. 
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B.  EQUALIZATION 


A  channel  accounts  for  the  effects  of  a  transmission  medium  as  a  signal  travels 
from  a  source  to  a  destination.  In  the  channel,  a  signal  can  be  amplified,  attenuated, 
dispersed  in  time,  and  distorted  as  a  function  of  the  characteristics  of  the  physical 
medium,  which  may  be  wireless  or  wired  along  with  miscellaneous  equipment.  As  a 
result,  every  channel  has  a  frequency  response  characteristic  that  will  affect  the  signal  as 
it  passes  and  will  distort  the  signal  from  its  original  form.  The  purpose  of  equalization  is 
to  nullity  the  effects  of  the  channel  characteristics  and  preserve  the  signal  as  close  as 
possible  to  its  original  form. 

Figure  7  is  a  schematic  model  of  a  channel  and  equalizer  as  used  in  a  variety  of 
applications.  In  this  work,  signal  acquisition  (sampling  and  binary  conversion)  is 
considered  part  of  the  channel.  The  discrete-time  signal  x[n]  is  the  observed  signal  at  the 
output  of  the  channel,  which  is  a  result  of  an  input  continuous-time  signal  s(t )  passing 
through  the  channel  and  sampled  at  the  receiver.  The  equalization  filter  attempts  to 
reconstruct  the  original  signal;  the  reconstructed  signal  s[n ]  is  desired  to  be  as  close  to 
s(nTs )  as  possible  where  Ts  is  the  sampling  interval. 


Figure  7.  Generic  Equalization  of  a  Signal,  s(t),  that  passes  through  a  Channel,  is  Digitized 
to  record  as  a  discrete -time  signal  x[n\  and  then  Equalized,  ■?[»] . 


If  a  channel’s  characteristics  are  known,  then  it  is  possible  to  effectively  remove 
the  impact  of  these  effects  by  utilizing  a  filter.  For  the  purpose  of  this  work,  all 
processing  was  perfonned  with  discrete-time  samples,  assuming  that  the  original  signal 
was  accurately  recorded  without  non-linearities,  represented  as  s[n],  and  discussed  further 
in  Chapter  III. 

The  channel  output  x\ri\  can  be  obtained  as  a  convolution  of  the  input  signal  s\ri\ 

and  the  channel  impulse  response  h[n\,  given  by 
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(2.16) 


x\n\  =  J'l  h[k]s[n  -  k]  n  =  0, 1, _ —  1 

k=0 

where  L  is  the  length  of  the  filter  representing  the  channel.  The  linear  operation  of 
convolution  can  be  accomplished  in  the  frequency  domain  as  a  product: 

X{co)  =  S{(o)H{(o)  (2.17) 

where  X(co)  is  the  Discrete  Time  Fourier  transform  of  the  received  signal  x\n\  and  I I(oj)  is 
the  channel  frequency  response.  Frequency  domain  implementation  of  convolution  using 
the  FFT  is  computationally  efficient  and  is  the  preferred  method  for  signal  analysis  in 
many  applications. 

The  system  model  for  equalization  of  a  signal  traversing  through  two  filters 
representing  the  channel  and  the  equalization  filter  is  shown  in  Figure  8.  The  input  signal 
to  the  channel  is  a  discrete-time  signal  s\n\.  In  this  model,  the  channel  is  assumed  to 
account  for  any  imperfections  in  the  data  acquisition  system.  Ideally,  the  equalization 
filter  will  have  a  frequency  response  that  is  the  inverse  of  the  channel.  In  other  words,  the 
overall  system  response  in  Figure  8  is  unity: 

H(a>)G(co)  =  1  where  G(co)  =  — —  (2.18) 

H(co) 

where  H{co)  is  the  channel  frequency  response  and  G(oj)  is  the  frequency  response  of  the 
equalization  filter. 


Figure  8.  A  Generic  Block  Diagram  of  Equalization  System  with  discrete  time  input  signal 
s\ri\.  The  Channel  and  Equalization  Filter  Blocks  are  represented  by  the  impulse 
responses  h[n]  and  g[/;]  respectively. 
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In  summary,  continuous-time  as  well  as  discrete -time  signals  can  be  processed  in 
either  the  time  domain  or  frequency  domain  through  the  utilization  of  various  transforms, 
such  as  the  FFT  and  DWT.  Although  the  Fourier  transform  is  considered  the  more 
traditional  method,  the  characteristics  and  versatility  of  wavelet  transforms  help  represent 
signals  in  multiple  ways.  Choices  in  wavelet  basis  function  and  the  number  of 
decomposition  levels  can  provide  a  collection  of  multi-band  waveforms  that  can  then  be 
filtered  or  equalized.  These  waveforms  can  also  be  processed  with  wavelet-specific 
techniques,  such  as  thresholding  to  de-noise  signals.  In  the  next  chapter,  these  techniques 
will  be  applied  to  the  NAWC-AD  aircraft  testing  process. 
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III.  SIGNAL  PROCESSING  OF  AIRCRAFT  TEST  WAVEFORMS 


This  chapter  outlines  the  overall  signal  processing  of  the  aircraft  testing  system 
and  highlights  the  characteristics  and  parameters  associated  with  the  measured  EMP 
waveforms.  It  will  briefly  describe  the  aircraft  test  process  and  detail  the  system  modeled 
to  be  used  for  analysis,  including  the  structure  of  the  parameters  available  for  the  system. 
This  provides  a  description  of  the  measured  EMP  waveforms  and  the  given  channel 
response  as  well  as  a  description  of  preprocessing  techniques  of  the  recorded  signal  and 
channel  response  infonnation  prior  to  equalization. 

A.  AIRCRAFT  TEST 

Considerable  work  has  been  done  to  ensure  that  military  aircraft  can  sustain 
various  debilitating  effects  that  could  occur  while  in  battle.  EMP  is  a  hazardous  and 
significant  phenomenon  because,  for  many  of  today’s  new  aircraft,  pilot  control  is  based 
on  electronic  signals  instead  of  traditional  hydraulic  control.  The  latter  is  impervious  to 
the  effects  of  electromagnetic  radiation  while  the  former  is  extremely  vulnerable  such 
that  an  EMP  could  disable  an  entire  aircraft. 

The  aircraft  EMP  testing  process  is  illustrated  in  Figure  9.  The  test  aircraft  is 
subjected  to  EMP  from  an  antenna  in  a  controlled  environment  [7].  The  test  waveform  in 
the  system  is  generated  by  an  EMP  antenna,  per  reference  MILSTD-464A,  and  then 
recorded  once  it  penetrates  the  aircraft.  Inside  the  aircraft,  there  are  multiple  sensors  that 
record  the  signal  at  designated  points.  These  measured  signals  are  subjected  to  the 
channel  effects  as  they  travel  through  the  recording  equipment  until  they  are  finally 
sampled  and  stored  in  digital  form. 
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channel  effects  of 
recording  equipment 


Figure  9.  Illustration  of  Aircraft  Testing  Process. 


Once  the  measured  EMP  waveform  has  been  captured  and  pre-processed,  it  is 
used  by  NAWC-AD  to  verify  that  the  signal  at  the  test  points  does  not  exceed  limits  of 
power  at  designated  frequencies.  By  knowing  the  characteristics  of  a  transmitted  test- 
EMP  waveform  and  the  waveform  acquired  inside  the  aircraft  at  designated  test  points, 
the  frequency  response  of  the  aircraft  structure,  as  a  medium,  can  be  detennined  and 
evaluated.  Military  aircraft  perfonnance  specifications  define  the  required  allowable 
levels  at  specific  frequencies;  therefore,  it  is  necessary  to  ensure  that  the  waveforms  at 
each  of  the  aircraft  test  points  are  accurately  acquired  for  processing.  This  detennines  the 
overall  effectiveness  of  the  aircraft  outer  shell  and  the  quality  of  radiation  hardening  in 
order  to  mitigate  the  damaging  effects  of  the  EMP  waves  [8]. 

Each  piece  of  equipment  in  the  acquisition  suite  has  non-ideal  frequency 
characteristics  and  other  limitations  that  will  distort  the  measured  waveforms. 
Consequently,  it  is  necessary  to  perfonn  equalization  of  the  imperfect  effects  of  the 
acquisition  hardware  to  recover  the  true  received  waveform  at  each  test-point.  This 
defines  the  objective  to  be  addressed  in  this  thesis:  to  ensure  that  the  recorded  EMP 
waveform  is  as  close  an  estimate  as  possible  of  the  waveform  that  would  have  been 
acquired  without  the  effects  of  the  acquisition  equipment,  and  to  consider  whether  the 
process  can  be  improved  with  other  techniques. 
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B. 


EMP  SIGNAL  COLLECTION  MODEL 


The  aircraft  testing  process  can  be  modeled  as  a  block  diagram  as  shown  in  Figure 
10.  Despite  the  complexity  of  the  aircraft  test  process,  the  essential  waveforms  required 
are  s(t),  x[n\,  and  s[n\ .  The  signal  s{t)  is  the  desired  signal,  x[n]  is  the  recorded  signal, 
and  .?[«]  is  the  estimate  of  the  desired  signal.  The  waveform  s(t)  is  considered  the  signal 
of  interest  since  it  has  penetrated  through  the  aircraft  chassis,  resultant  of  a  controlled 
transmitted  wideband  EMP  wave.  The  information  contained  in  the  waveform  provides 
the  analysis  for  the  aircraft  performance.  This  signal  then  traverses  through  various 
pieces  of  recording  and  transfer  equipment  until  it  is  ultimately  digitized  and  stored.  This 
discrete-time  recorded  waveform  is  designated  as  x[n ].  Finally,  the  signal  is  processed 
using  equalization  techniques  to  remove  the  effects  of  the  acquisition  equipment  and 
represented  as  the  signal  estimate .?[/?] . 


Channel 

Figure  10.  Block  Diagram  of  the  Aircraft  Test  Model. 


Although  there  are  multiple  pieces  of  equipment  that  comprise  the  acquisition 
suite  and  recording  equipment,  the  overall  system  from  sensors  to  storage  can  be 
expressed  as  a  consolidated  single  channel. 

C.  EQUALIZATION  OF  THE  ACQUIRED  WAVEFORM 

Although  the  test-EMP  waveforms  and  aircraft  chassis  characteristics  are  the 
factors  that  make  up  the  basis  of  NAWC-AD’s  testing,  they  are  irrelevant  to  the 
equalization  and  signal  processing  problem  discussed.  Figure  1 1  models  the  process  from 
the  desired  signal  to  the  equalized,  approximated  signal. 
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Due  to  the  nature  of  the  comparison,  it  is  not  necessary  to  represent  s(t )  as  a  continuous 
waveform  for  the  system.  The  equipment  and  the  signal  processing  techniques  are  digital 
and  discrete,  and  it  is  assumed  that  the  sampling  losses  of  s(t)are  built  into  the  NAWC- 
AD  evaluation  criteria.  As  a  result,  s(t)  in  this  system  is  best  represented  as  a  sampled 
waveform  .v[/?]  with  sampling  time  Ts  where  s\n\  =  s(nTs ) . 


S[«]  *[«]  s[ft] 


Figure  1 1 .  Equalization  Model. 


The  channel  and  the  equalization  filter  can  both  be  represented  as  systems  with 
response  infonnation  as  discussed  in  Chapter  II.  The  channel  characteristics  are  given 
quantities  based  on  equipment  specification  and  test  results.  The  channel  information  is 
provided  by  NAWC-AD  in  the  frequency-domain  and  will  be  as  denoted  as  H(co).  From 
this  a  priori  knowledge  of  the  channel,  it  is  possible  to  construct  an  equalization  filter, 
G(co)  to  counteract  the  effects  of  the  recording  equipment.  Applying  the  discussion  from 
Chapter  II,  it  can  be  shown  that 

X(co)  =  S(cq)H(co) 

(3.1 

S(co)  =  X{co)G{(o )  =  S(co)H(co)G(co) 


Consequently,  in  order  to  ideally  equalize  the  system  such  thatkfin)  =  S(co) ,  we 

desire  G(a> )  =  — - — .  Using  this  relationship  and  the  provided  channel  response  H(co) 
H(co) 

from  NAWC-AD,  it  was  possible  to  derive  a  G(co)  for  experimentation.  In  discrete -time, 
this  vector  must  have  unifonn  frequency  distribution  as  well  as  be  of  the  same  length  as 
the  signal  vector.  Both  these  conditions  were  required  to  be  met  through  pre-processing 
due  to  the  fonnat  of  the  H(co)  information  discussed  later. 
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D.  MEASURED  SIGNAL  CHARACTERISTICS 


Unfortunately,  it  is  not  possible  to  gather  information  on  s(t)  because  it  is  the  true, 
unsampled,  and  unrecorded  signal.  Only  x\n\  can  be  made  available  since  it  is  the 
recorded  signal  at  the  aircraft  test  points;  therefore,  it  is  necessary  to  use  a  representative 
signal  for  the  purposes  of  signal  analysis. 

The  given  acquired  signal,  x\n\,  is  sampled  with  a  fixed  sampling  rate  in  the  range 
of  250  MHz  to  1 .5  GHz,  dependent  upon  the  equipment  used.  The  waveform  is  a  real¬ 
valued  signal  and  varies  in  length  between  9,000  and  25,000  points  per  vector.  As  a 
result,  s\n\  and  v[/;]  will  have  the  same  discrete  characteristics  as  x\n\.  Visually,  the 
signal  reflects  the  shape  of  a  transient  waveform  as  seen  in  Figure  12  (a)  though  (d).  The 
reader  may  note  that  the  waveform  shape  is  not  consistent  as  is  a  function  of  the  location 
of  test  points  within  the  aircraft. 


Time  in  sec  x10"5 
(a) 


Time  in  sec 

(b) 


Time  in  sec 

(c) 


x  10 


Time  in  sec 

(d) 


xIO 


Figure  12.  Example  Recorded  Waveforms  at  Various  Aircraft  Test  Points,  (a) 

Wavefonn  1.  (b)  Wavefonn  2.  (c)  Waveform  3.  (d)  Waveform  4. 
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There  appears  to  be  a  low  level  of  additive  noise  in  the  recorded  signal  x\ri\.  It  is 
not  evident  whether  if  this  noise  is  equipment  noise  or  background  noise  resident  during 
the  radiating  of  the  aircraft.  Accordingly,  the  recorded  signal  can  be  represented  as  either 

x[n\  =  s\ri\  *  h[n]  +  u[/;] 


or 


x[n]  =  (^[//]  +  w\n\)  *  /?[/?] . 

This  noise  is  present  in  the  recorded  signal  but  may  not  be  part  of  the  signal  ,s'[n].  It  is 
assumed  that  the  noise  is  undesirable  and  not  part  of  the  ideal  system;  therefore,  the 
results  could  be  more  accurate  if  it  were  not  present. 


E.  MEASURED  CHANNEL  RESPONSE 

Provided  by  NAWC-AD  were  frequency  responses  for  the  acquisition  equipment 
used  in  collecting  the  EMP  signals.  The  length  of  the  equipment’s  channel  vectors  varied 
from  800-1,200  points  with  frequency  values  spanning  from  200  MHz  up  to  1  GHz.  The 
channel  response  is  complex-valued,  and  the  frequencies  do  not  have  unifonn  frequency 
spacing  conducive  to  Fourier  processing.  Each  piece  of  equipment  was  sampled  using 
multiple  sampling  frequencies  and  each  channel  was  sampled  in  different  sampling 
patterns. 

In  Figure  13,  (a)  and  (c)  illustrate  two  different  channel  spectra  provided  by 
NAWC-AD.  The  spectrum  in  Figure  13(a),  “8694. cal,”  is  a  vector  of  1 196  magnitude 
values  across  a  range  of  1  GHz.  The  minimum  spacing  between  these  frequencies  was  28 
Hz,  with  the  maximum  difference  being  4.166  MHz.  Figure  13(b)  illustrates  the 
irregularity  of  sampling  by  plotting  the  difference  in  frequencies  between  each  successive 
frequency  values  across  the  entire  H(co)  vector.  It  appears  that  there  is  a  step  like 
sampling  pattern  based  upon  the  location  in  the  frequency  spectrum.  Similarly,  Figure 
13(c)  illustrates  the  H(co)  characteristics  for  the  “wbal.cal”  vector  which  is  made  up  of 
1 196  frequency  values.  It  has  differences  in  frequency  spacing  ranging  from  5  Hz  to  2.79 
MHz,  although  its  behavior  is  exponential  in  shape  as  seen  in  Figure  13(d). 
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Magnitude  Magnitude 


As  discussed  in  Chapter  II,  discrete  channel  responses  have  equal  and  uniform 
spacing  throughout  the  spectral  range  of  interest  such  that  there  can  be  a  one-to-one 
correspondence  with  the  time  domain.  Also,  the  given  channel  frequency  response 
represents  only  the  positive  frequencies  from  0  <  co  <  n  .  Lastly,  the  measured  channel 
response  provided  is  sparsely  defined  where  800-1,200  discrete  points  represent  an  entire 
1GHz  bandwidth. 


Frequency  (MFIz)  Frequency  Index 


(c)  (d) 

Figure  13.  Channel  Responses  and  the  Sampling  Rates  of  Channel,  (a)  Channel 

Frequency  Response  of  “8694. cal.”  (b)  Difference  in  Frequency  Values  of 
Successive  Frequency  Index  Points  in  Channel  File  “8694. cal.”  (c)  Channel 
Frequency  Response  of  “wbal.cal.”  (d)  Difference  in  Frequency  Values  of 
Successive  Frequency  Index  Points  in  Channel  File  “wbal.cal.” 
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As  a  result,  there  is  required  processing  of  the  equipment  frequency  response 
vectors  to  conform  to  uniform  frequency  spacing  of  the  discrete  Fourier  domain. 
Additionally,  the  spectrum  will  require  having  the  digital  frequency  values 
from  —it  <  <x>  <  0  .  As  a  note,  irregular  sampling  rates  can  be  advantageous  if  gathered 
appropriately  for  a  few  reasons.  Since  Fourier  transforms  require  uniform  sampling  rates 
throughout  the  spectrum,  high  sampling  frequencies  and  wide  bandwidths  can  take  up 
large  amounts  of  storage  as  well  as  require  significant  amounts  of  processing  time.  There 
is  the  potential  to  exploit  these  multi-rate  characteristics  in  a  more  optimal  manner  to 
support  using  wavelet  techniques  versus  Fourier  transforms. 

F.  SIGNAL  PRE-PROCESSING 

There  are  three  types  of  pre-processing  that  NAWC-AD  performed  prior  to 
equalization.  These  are  signal  extraction,  DC  bias  removal,  and  linear  interpolation  of 
channel  characteristics.  Before  it  can  be  determined  whether  more  effective  methods  can 
be  employed,  existing  processes  must  be  discussed. 

1.  Signal  Extraction 

The  technique  of  signal  extraction  in  the  time  domain  is  utilized  such  that  signal 
processing  only  occurs  on  information  where  the  signal  is  present.  There  is  the  potential 
for  adversely  impacting  the  information  by  utilizing  a  portion  of  data  in  which  only  noise 
is  present.  Accordingly,  it  is  beneficial  to  eliminate  recorded  information  that  doesn’t 
contain  the  signal  both  before  and  after  an  EMP  signal  pulse.  Signal  extraction  is  done  by 
NAWC-AD  in  the  time  domain  through  manual  processing  and  through  visual  inspection 
of  the  recorded  waveform.  In  each  case,  the  information  preceding  the  signal  impulse  is 
manually  removed  by  the  operator,  and  the  data  at  the  end  of  the  recorded  waveform  are 
removed  after  it  appears  as  though  no  signal  is  present  in  the  recorded  information.  This 
method  eliminates  consistency  from  the  test  and  evaluation  process  because  it  is  manual 
and  relies  on  a  user  interface  decision.  In  this  manner,  often  the  same  data  can  produce 
varying  degrees  of  processing  effectiveness  due  to  the  proficiency  of  the  operator. 

For  a  signal  with  no  noise,  it  will  have  no  average  energy  where  the  signal  is  not 
present  and  a  positive  value  of  energy  wherever  the  signal  is  present.  Any  additive  white 
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noise  theoretically  will  create  a  noise  floor.  The  average  energy  for  a  signal  over  a  given 
data  length  interval  can  be  found  as  shown 

EAvM  =  \-j^\x[n-if  (3.2) 

N  i=N- 1 

where  N  is  length  of  the  signal  vector.  In  portions  where  the  signal  is  present,  the  average 
energy  of  the  signal  is  additive  to  the  noise  floor. 

In  order  to  perfonn  this  in  an  automated  manner,  it  was  necessary  to  determine  a 
threshold  value  that  detennines  where  the  signal  begins  and  where  it  ends.  The  threshold 
is  calculated  based  on  the  average  energy  found  at  the  beginning  and  end  of  the  measured 
signal.  For  this  algorithm  to  be  successful,  it  must  be  assumed  and  required  that  there  is 
no  signal  at  the  beginning  or  end  of  the  measured  signal. 

Figure  14  (a)  shows  a  sample  waveform  that  has  been  recorded,  that  has  sample 
values  at  the  start  and  end  of  the  information  that  do  not  contain  the  desired  signal.  When 
the  average  energy  is  calculated  and  plotted,  as  shown  in  14(b),  the  noise  floor  can  be 
noted  at  the  start  and  end  of  the  signal  energy,  and  the  location  of  the  actual  signal  is 
evident.  When  extraction  techniques  are  applied  to  the  average  energy  plot,  as  shown  in 
14(c),  by  using  a  threshold  value,  the  portions  of  the  signal  assumed  to  have  no  signal 
elements  are  removed,  leaving  only  the  length  of  the  infonnation  containing  the  signal. 
Lastly,  14(d)  shows  the  extracted  signal  with  non-signal  portions  removed. 
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Figure  14.  Example  of  Signal  Extraction,  (a)  Original  Waveform,  (b)  Time  Average 

Energy  of  Waveform,  (c)  Range  of  Time  Average  Energy  After  Signal  Points 
Below  Threshold  Have  Been  Removed,  (d)  Final  Extracted  Signal. 


For  this  process,  if  the  threshold  is  too  large  in  magnitude,  a  portion  of  the  signal 
could  be  lost  which  may  have  an  adverse  affect  on  the  processing.  If  the  threshold  is  too 
small,  it  is  possible  that  the  noise  is  not  completely  removed  and  a  portion  of  the  recorded 
samples  that  have  no  presence  of  the  signal  is  retained.  This  additional  noise  could  distort 
results.  The  selection  of  the  appropriate  threshold  value  may  have  to  be  done  by  trial  and 
error. 

Figure  15  shows  the  effects  of  using  different  threshold  values.  Figure  15(a)  is  the 
actual  measured  wavefonn  using  no  threshold  value.  It  characterizes  the  case  where 
portions  of  the  recorded  waveform  containing  no  energy  of  the  signal  remain.  With  a 
threshold  value  of  1,  as  in  Figure  15(b),  the  extracted  signals  appear  to  have  the  correct 
balance  of  signal  removal  that  is  expected.  In  Figures  15  (c)  and  (d),  the  threshold  value 

was  selected  to  be  too  large,  and  significant  portions  of  the  signal  are  removed.  By  visual 
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Amplitude  Amplitude  Amplitude  Amplitude 


inspection,  a  threshold  value  close  to  1  appears  to  be  optimum  for  this  wavefonn  as  it  is 
extracting  the  entire  signal  present  in  the  data  with  the  potions  of  the  noise  removed. 


(d) 


Figure  15.  Signal  with  Varying  Levels  of  Energy  Threshold  Values,  (a)  Threshold  = 

0.1.  (b)Threshold  =  1.  (c)  Threshold  =  20.  (d)  Threshold  =  80. 


2.  DC  Bias  Cancellation 

System  electrical  components  may  contain  a  DC  bias  in  their  equipment  that 
would  shift  the  values  of  the  entire  waveform.  This  is  an  undesired  effect  of  the 
acquisition  electronic  equipment  and  may  be  removed.  The  technique  of  averaging  is 
performed  to  eliminate  the  effects  of  the  DC  bias  in  system  equipment.  In  an  enviromnent 
where  no  signal  is  present,  it  is  assumed  that  the  average  value  of  energy  would  be  zero 
as  it  is  with  white  noise.  When  signals  travel  though  electronic  equipment,  an  artificial 
DC  gain  can  manifest  itself  onto  the  signal.  DC  averaging  is  the  technique  to  determine 
that  value  such  that  it  can  be  removed  from  all  the  data.  The  current  method  that  NAWC- 
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AD  uses  for  DC  averaging  is  to  take  the  first  20-150  waveform  values  of  the  recorded 
signal  and  compute  the  average  value  of  these  points.  This  value  is  then  removed  from  all 
signal  sample  values. 

As  far  as  the  experiments  are  concerned,  it  was  not  possible  to  verify  whether  the 
cancellation  of  DC  bias  actually  improved  performance.  Because  the  DC  bias  is 
introduced  into  the  channel  in  the  acquisition  suite,  the  noise  contribution  would  need  to 
be  artificially  added,  only  to  be  removed.  With  no  true  signal  s[n\  available,  it  is  not 
possible  to  verify  this  improvement.  Therefore,  there  are  no  results  provided  for  this 
topic. 


3.  Linear  Interpolation  of  the  Channel  Frequency  Response 

As  mentioned,  the  channel  response  is  provided  in  non-uniform  sampling 
intervals.  In  order  to  provide  equalization  of  the  discrete  waveforms  in  the  method 
discussed  above,  both  the  signal  vector  and  the  channel  must  have  an  equal  number  of 
points  and  the  same  digital  frequency  in  order  to  perform  processing.  In  the  channel 
information^  raw  form,  there  is  a  disparity  between  the  length  of  the  signal  vectors  and 
the  channel  frequency  response.  Therefore,  the  channel  response  requires  preprocessing 
before  it  can  be  applied  to  the  signal. 

It  is  necessary  to  fill  the  channel  spectrum  with  magnitude  values  for  each  of  the 
data  frequency  intervals  or  bins.  This  is  accomplished  through  linear  interpolation  of  the 
channel  response.  A  uniformly  sampled  data  set  of  length  N,  sampled  at  a  rate  of  Fs,  will 
have  an  FFT  of  length  N  with  frequency  bins  spaced  FJN  Hz  apart.  The  channel 
responses  for  the  given  equipment  do  not  have  such  a  uniform  sampling  as  described 
above. 

In  order  to  equalize  the  measured  waveform,  a  magnitude  value  must  be  found  for 
the  channel  response  at  each  discrete  frequency  of  the  data.  As  shown  in  Figure  16,  a 
magnitude  is  linearly  interpolated  from  successive  points  in  the  channel  frequency 
response  for  frequency  values  dictated  by  the  measured  EMP  waveform  frequency  bins. 
This  process  is  iterated  for  each  FFT  frequency  bin  so  that  the  resultant  channel  response 
will  have  an  equal  number  of  points  with  corresponding  magnitude  values  at  each 
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frequency  represented  in  the  signal’s  frequency  response.  Those  values  of  the  channel 
response  that  extend  past  half  of  the  sampling  frequency  of  the  data  are  disregarded. 


A 


channel 


A 


data 


ft 


channel+1 


Frequency 


Figure  16.  Linear  Interpolation  of  Channel  Response.  Frequencies  f channel  and 

fchannei+i  Represent  Successive  Points  in  the  Channel  Spectrum.  Frequency  fdata 
Represents  Location  of  the  Frequency  Bin  of  Data  in  the  Frequency  Domain. 

This  chapter  discussed  the  aircraft  testing  scheme  as  performed  by  NAWC-AD 
for  EMP  pulses  as  well  as  discussed  the  model  for  which  they  equalize  the  raw  results.  It 
introduced  the  concept  of  equalization  as  applied  to  the  NAWC-AD  process  of  removing 
acquisition  effects  to  the  system.  This  chapter  also  presented  the  measured  signal 
characteristics  and  the  measured  channel  responses,  noting  challenges  and  preprocessing 
requirements  prior  to  perfonning  signal  equalization.  Lastly,  this  chapter  described  three 
pre-processing  techniques  that  NAWC-AD  implemented  prior  to  equalizing  the  signal 
data.  The  next  chapter  will  present  four  equalization  techniques  using  Fourier  and 
wavelet  transfonns,  as  well  as  the  implementation  of  de-noising  and  thresholding. 
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IV.  EQUALIZATION  METHODS  AND  WAVELET  SCHEMES 


The  previous  chapters  discussed  the  system  and  parameters  pertinent  to  measured 
EMP  waveforms  and  the  channel  as  well  as  pre-processing  and  equalization  techniques 
that  may  be  applied  to  processing  the  measured  information.  Three  different  equalization 
methods  were  designed  to  process  the  recorded  EMP-signals  using  wavelet  transfonns  to 
acquire  results.  The  EMP  signals  were  also  processed  using  the  Fourier  transform 
method  currently  used  by  NAWC-AD  for  comparison.  This  chapter  details  the  algorithms 
and  steps  taken  to  equalize  the  waveforms  and  to  reconstruct  the  signal  received  at  the 
aircraft  test-points. 

A.  EQUALIZATION 

Using  the  principles  described  in  Chapter  II  and  the  knowledge  of  the  system,  a 
system  model  was  required  to  simulate  the  channel  and  filters.  Using  recorded  data  from 
NAWC-AD,  it  was  possible  to  implement  a  model  of  the  aircraft  testing  and  recording 
system  in  MATLAB  and  apply  wavelet-based  signal  processing  techniques  in  place  of 
those  currently  used  by  NAWC-AD.  The  remainder  of  this  chapter  summarizes  the 
techniques  employed  for  this  analysis. 

1.  Equalization  using  Wavelets 

To  implement  equalization  using  wavelet  transforms  and  perform  filtering  at  the 
sub-band  level,  a  generic  model  of  the  system  was  required  to  employ  wavelet  techniques 
as  shown  in  Figure  17.  The  process  requires  wavelet  decomposition,  followed  by 
filtering,  and  then  synthesis  of  the  filtered  wavefonn. 


s[«] 


Figure  17.  System  model  using  Wavelet  Processing  for  Equalization. 
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Focusing  on  the  equalization  process,  a  filter  bank  structure  as  shown  in  Figure  18 
illustrates  the  equalization  process.  The  input  signal  x[n\  is  decomposed  into  sub-band 
components  using  the  filters  h'  and  g'  discussed  in  Chapter  II,  filtered  for  equalization, 
and  then  reconstructed  into .?[/?] .  Each  individual  filter  lengths  of  equalization  filters  1-4 
are  matched  to  the  length  of  sub-band  signals  for  processing.  All  filters  were  developed 
and  implemented  using  MATLAB. 


.v[n] 


Figure  18.  Structure  of  Wavelet  Analysis  and  Synthesis  with  Equalization  Filtering. 


The  measure  of  system  performance  for  the  purpose  of  comparing  the  algorithms 
developed  in  this  chapter,  the  Mean  Squared  Error  (MSE)  between  the  desired 
waveforms  s[n]  and  the  estimated  waveform s[n\,  is  computed  as  [9] 

1  N- 1  2 

^ms=— 2>KMKI)  (41) 

n= o 

where  N  is  the  numbers  of  signal  values,  and  the  desired  waveform,  s[n\,  is  compared 
against  the  estimated  signal,  s[n\ .  Mean  squared  error  is  one  of  many  possible  measures 
of  comparison,  but  it  is  the  only  one  utilized  to  compare  the  results  of  the  experiments 
reported  here. 

2.  Equalization  Techniques 

Three  algorithms  were  derived  to  accomplish  this  equalization  processing.  All 
techniques  were  based  on  wavelet  decomposition  of  recorded  EMP  waveforms,  but 
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varied  in  types  of  wavelet  functions,  filter  orders,  and  preprocessing  techniques.  The 
three  best  perfonning  methods,  as  well  as  the  technique  used  by  NAWC-AD,  were 
evaluated  and  compared.  The  fourth  method  is  the  model  of  equalization  employed  by 
NAWC-AD. 


a.  Method  1  -  Time  Domain  Filtering  of  Sub-band  Time  Domain 
Signals  with  Inverse  HR  Filter 

The  first  method  of  equalization  is  an  algorithm  that  decomposes  the  EMP 
signal  x[n]  into  sub-band  time  domain  signals,  and  filters  each  sub-band  signal  by  an  HR 

filter  constructed  using  the  Yule-Walker  method  having  the  fonn — — . 

B(z) 

As  shown  in  Figure  19,  the  recorded  EMP  signal,  x[n],  is  decomposed 
using  the  DWT  into  approximation  and  detail  coefficients.  Each  set  of  coefficients  is  then 
selectively  reconstructed  as  discussed  in  Chapter  II  into  time  domain  sub-band  signals  of 
appropriate  length  and  sampling  frequency  of  the  data.  These  multiple  signals  are  then 
filtered  and  equalized  in  the  time  domain.  The  resultant  time-domain  sub-band  signals  are 
decomposed  back  into  their  respective  coefficients,  and  then  reconstructed  back  into  a 
full  spectrum  time  domain  signal.  During  decomposition,  sub-band  coefficients  outside 
the  filtered  sub-bands  are  discarded  and  coefficients  are  set  to  zero. 


Figure  19.  Method  1:  Processing  of  Signal  Waveforms  Using  Wavelet-Based 

Techniques. 
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b.  Method  2  -  Frequency  Domain  Filtering  using  Channel 
Frequency  Response 

The  second  method  of  equalization  occurs  in  the  frequency  domain  versus 
the  time  domain.  Like  Method  1  (see  Figure  20),  the  algorithm  decomposes  the  signal 
into  coefficients  and  reconstructs  the  sub-band  signals  into  the  time  domain.  It  then  takes 
the  Fourier  transform  of  each  sub-band  signal  and  then  filters  each  sub-band  signal  by 
dividing  the  frequency  response  of  the  channel.  The  resultant  frequency  domain  signals 
are  transformed  back  to  the  time  domain  using  the  inverse  Fourier  transform,  where  the 
time-domain  sub-band  signals  are  decomposed  into  wavelet  coefficients.  The  respective 
coefficients  are  then  used  to  reconstruct  the  complete  signal  back  to  the  time  domain. 


Figure  20.  Method  2:  Processing  of  Signal  Waveforms  Using  Wavelet-Based 

Techniques,  and  Fourier  Transforms. 


c.  Method  3  -  Time  Domain  Filtering  of  Wavelet  Coefficients  using 
Inverse  HR  Filter 


The  third  algorithm,  as  shown  in  Figure  21,  decomposes  the  EMP  signal 
x[n]  into  wavelet  approximation  and  detail  coefficients.  Then  each  of  the  sub-band  filter 


coefficient  sets  is  filtered  by  the  HR  equalization  filter 


Az) 

B(z) 


designed  based  on  the 


channel  response  infonnation.  The  signal  is  then  reconstructed  using  the  sub-band 
coefficients.  The  difference  between  this  procedure  and  that  of  Method  1  is  that  the 
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coefficients  are  not  selectively  reconstructed  back  to  their  time  domain  forms.  The 
signal’s  sub-band  components  of  variable  lengths  are  fdtered  accordingly. 


Figure  2 1 .  Method  3:  Processing  of  Signal  Waveforms  Using  Wavelet-Based 

Techniques  and  HR  Filtering. 


d.  Method  4 -NA  WC-AD  Method 

The  NAWC-AD  method  filters  the  recorded  signal,  x\n\,  in  the  time 
domain  by  using  the  vector  multiplication  of  the  channel’s  inverse  response  as  shown  in 
Figure  22.  The  time  domain  signal  is  then  recovered  by  using  the  IFFT.  This  method  was 
discussed  in  Chapter  II. 


Figure  22.  Method  4:  Processing  of  Signal  Waveforms  Using  Fourier  Based 

Transforms  as  Performed  by  NAWC-AD. 
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B.  SUB-BAND  FILTERING  AND  NOISE  REMOVAL 

Although  de-noising  is  a  powerful  technique  available  using  wavelet  transforms,  a 
detailed  evaluation  was  not  explored  in  this  thesis.  However,  as  a  demonstration  of  de- 
noising  techniques,  a  soft  thresholding  de-noising  filter  using  a  threshold  based  on  Stein’s 
unbiased  risk  calculations  (provided  by  MATLAB)  was  applied  at  each  sub-band  to 
produce  results  for  comparison.  As  shown  in  Figure  23,  the  process  for  de-noising  occurs 
after  signal  decomposition  and  before  equalization.  The  threshold  calculation  noted 
above  is  standard  across  the  entire  signal. 


Figure  23.  De-noising  Model  for  Processing  EMP  Wavefonns. 


Further  investigation  would  be  necessary  to  detennine  the  best  method  for  each 
given  channel,  signal  type,  and  wavelet  basis  function. 

1.  System  Model  with  Additive  White  Gaussian  Noise 

To  further  analyze  the  proposed  equalization  techniques,  an  investigation  was 
conducted  into  the  performance  of  these  Fourier  and  wavelet  transform  methods  when 
there  is  noise  present.  A  system  model  with  additive  white  Gaussian  noise  of  varying 
levels  added  to  the  signal  was  developed  as  shown  in  Figure  24.  Noise  is  introduced  into 
the  system  by  adding  it  to  the  recorded  signal,  x[n\. 
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w[n\ 

Figure  24.  System  Diagram  with  the  Addition  of  White  Noise,  w\n\. 


Using  different  methods  of  equalization  described  above,  the  resultant  waveform 
s[n\  was  compared  to  the  original  noise-free  waveform  over  a  wide  range  of  w[n]  power 
levels.  This  further  expands  the  analysis  of  this  experiment  by  comparing  the 
performance  of  the  algorithms  in  noise-free  and  noisy  environments  and  will  be  further 
discussed  as  part  of  the  results  in  Chapter  V. 

This  chapter  provided  the  specific  models  that  were  implemented  to  verify 
whether  or  not  wavelet-based  methods  could  outperform  existing  methods.  It  outlined 
methods  and  configurations  used  for  equalization.  Finally,  it  addressed  wavelet 
thresholding  and  the  potential  for  data  de-noising.  The  next  chapter  documents  the  results 
of  the  system  simulations  and  experimental  results. 
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V.  SIMULATION  RESULTS 


The  previous  chapters  discussed  many  of  the  signal  processing  principles  and 
techniques  that  apply  to  the  NAWC-AD  EMP  waveform  processing.  An  experiment  was 
defined  in  order  to  test  proposed  alternative  techniques.  The  experiment  was  conducted 
with  five  types  of  wavelets  using  the  four  methods  discussed  in  Chapter  IV.  Comparison 
between  three  and  five  levels  of  wavelet  decomposition  was  used  since  early 
experimental  analysis  showed  that  perfonnance  degraded  as  the  number  of  levels 
exceeded  five.  Additionally,  each  set  of  data  was  processed  using  signal  extraction  and 
de-noising  techniques  as  well.  This  chapter  presents  results  that  characterize 
performance  of  the  various  techniques  compared  to  the  present  method  of  signal 
processing.  As  discussed  in  Chapter  IV,  the  measure  of  performance  is  the  mean  squared 
error  between  the  input  data  and  the  equalized  estimate  of  the  measured  EMP  signal. 

Since  the  FFT-based  NAWC-AD  technique,  Method  4,  is  the  currently  used  method  for 
processing  the  EMP  signals,  all  proposed  techniques  will  be  measured  against  it. 

There  were  two  primary  experiments  conducted  for  this  thesis.  The  first  involved 
a  comparison  of  the  different  wavelet-based  techniques  for  equalization  of  the  system 
with  the  given  signal  information  provided  by  NAWC-AD.  The  second  experiment 
investigated  the  performance  of  these  methods  with  the  addition  of  additive  noise. 

A.  SIMULATION  MODEL  FOR  ANALYSIS  OF  SYSTEM  WITH  NO 

ADDITIVE  NOISE 

There  were  two  models  used  for  analysis  in  the  experiment.  All  models  employed 
an  HR  filter  representing  the  acquisition  equipment  channel  constructed  from  the  channel 
frequency  response  information  provided  by  NAWC-AD.  The  filter  is  implemented  in 
time  domain  to  keep  uniformity  in  its  implementation.  Figure  25  depicts  the  basic  process 
of  the  experiment. 
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Figure  25.  Basic  Model  of  System  Experiment. 


HR  Filter 


In  addition,  as  described  in  Chapter  III,  the  technique  of  signal  extraction  was 
noted  as  a  mechanism  of  improving  perfonnance  and  is  incorporated  into  the  model  as 
can  be  seen  in  Figure  26. 


Figure  26.  Basic  Model  of  System  Implementing  Signal  Extraction  Technique. 

Each  of  these  models  supports  analysis  of  different  types  of  wavelets  processing 
techniques  against  the  Fourier  transform  based  processing,  as  well  as  performance 
comparison  of  different  wavelets  against  each  other.  The  signal  used  for  all  experiments, 
with  exception  to  a  waveform  to  waveform  comparison  discussed  below,  used  Waveform 
1  depicted  in  Figure  12(a).  The  types  of  wavelets  used  in  the  process  were  arbitrarily 
selected  and  are  Haar,  Meyer,  Coiflet,  Daubechie’s  Order  3,  and  Daubechie’s  Order  15 
filters.  In  this  experiment,  3-5  levels  of  wavelet  decomposition  were  evaluated  as  it  was 
noted  during  experimentation  that  they  provided  the  best  results.  Mean  Squared  Error 
(MSE)  as  defined  in  Chapter  IV  was  the  basis  of  comparison.  The  technique  of  de- 
noising  is  compared  in  each  model;  de-noising  occurs  as  a  function  after  the  signals  have 
been  decomposed  into  wavelet  sub-bands.  A  comparison  of  these  results  will  follow  for 
each  model. 
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1.  Equalization  with  No  Signal  Extraction 

The  first  set  of  results  is  the  comparisons  of  the  mean  squared  error  values  with 
no  de-noising  as  shown  in  Figure  25.  It  is  a  direct  comparison  of  wavelet  methods, 
Methods  1-3,  against  the  Fourier  transform  technique,  Method  4,  using  different  wavelet 
basis  functions  and  various  levels  of  decomposition.  As  noted,  Waveform  1,  a  wideband 
waveform  of  16,384  points,  was  used  following  the  process  model  of  Figure  25  without 
signal  extraction.  Each  of  the  four  methods  of  equalization  was  discussed  in  Chapter  IV, 
and  wavelet  decomposition  is  provided  to  the  3rd,  4th  and  5th  levels.  The  channel 
frequency  response  used  in  these  models  was  based  on  “8694. cal”  depicted  in  Figure 
13a.  It  is  a  1,196-point  Fourier  domain  response,  non-linearly  spaced,  from  0-1  GHz.  This 
response,  therefore,  required  linear  interpolation  of  its  frequency  values  as  discussed  in 
Chapter  III,  Section  F.  Additionally,  the  mean  squared  error  values  are  expressed  in 
decibels  (dB)  since  the  absolute  values  are  relatively  small.  Results  of  these  simulations 
can  be  seen  in  Table  1. 
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Haar  Wavelet 


Level 

Method  1 

Method  2 

Method  3 

Method  4 

3 

-89.62654218 

-64.18408051 

-79.94143096 

-129.38375380 

4 

-88.02121882 

-64.18403540 

-76.75129362 

-129.38375380 

5 

-86.68836043 

-64.18404263 

-74.06419920 

-129.38375380 

Daubechies  Order  3  Wavelet 


Level 

Method  1 

Method  2 

Method  3 

Method  4 

3 

-89.14776518 

-64.17712981 

-79.74522967 

-129.38375380 

4 

-87.93476348 

-64.17463116 

-76.59660948 

-129.38375380 

5 

-86.36230323 

-64.17081395 

-73.82303546 

-129.38375380 

Daubechies  Order  15  Wavelet 


Level 

Method  1 

Method  2 

Method  3 

Method  4 

3 

-92.69689015 

-64.18181895 

-79.74779153 

-129.38375380 

4 

-91.12815898 

-64.17953860 

-76.54029373 

-129.38375380 

5 

-88.01481085 

-64.17318187 

-73.75515026 

-129.38375380 

Coiflet  Order  1  Wavelet 


Level 

Method  1 

Method  2 

Method  3 

Method  4  ; 

3 

-89.99558145 

-64.18042512 

-79.89040809 

-129.38375380 

4 

-88.63392335 

-64.17931756 

-76.68335041 

-129.38375380 

5 

-86.39083615 

-64.17317788 

-73.89602337 

-129.38375380 

Meyer  Wavelet 


Level 

Method  1 

Method  2 

Method  3 

Method  4 

3 

-78.08774551 

-63.52205199 

-76.54194364 

-129.38375380 

4 

-75.06072253 

-63.03518563 

-73.11895494 

-129.38375380 

5 

-72.76247718 

-62.46281427 

-70.42286280 

-129.38375380 

Table  1 .  Mean  Squared  Error  in  dB  Using  All  Three  Techniques  for  Various  Wavelet 
Basis  Functions  and  Multiple  Level  Decompositions  and  the  Fourier  method. 

Signal  is  Not  Extracted  from  Data. 

Overall,  the  perfonnance  of  Method  4  exceeds  that  of  all  other  methods.  The 
closest  performing  method  to  Method  4  is  Method  1  using  Daubechie’s  order  15  wavelet, 
with  three  levels  of  decomposition.  As  can  be  seen,  the  wavelet  basis  function  has  an 
impact  on  performance  as  well  as  the  level  of  decomposition.  The  performance  of 
Method  1  is  the  best  in  the  no  noise  comparison,  followed  by  Method  3  and  then  Method 
2.  Since  the  Method  1  filters  are  derived  directly  from  the  channel  frequency  response, 
this  is  expected.  Also,  the  performance  of  the  wavelet  methods  degrades  as  the  number 
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of  levels  is  increased  for  each  of  the  wavelet  families.  As  wavelet  basis  functions, 
Daubechie’s  order  15  wavelet  has  the  best  perfonnance  in  Method  1,  but  for  Methods  2 
and  3  the  Haar  wavelet  performs  the  best. 

2.  Equalization  with  Signal  Extraction 

The  next  set  of  results  captures  the  effectiveness  of  the  signal  extraction  algorithm 
described  in  Chapter  III,  Section  F.  For  consistency,  Waveform  1  and  channel  response 
“8694. cal”  were  used  again  with  the  same  preprocessing  and  equalization  methods  as  in 
the  above  discussion.  Table  2  captures  the  results  of  this  comparison  (same  results  as  in 
Table  1,  but  incorporating  signal  extraction  in  preprocessing  as  depicted  in  the  model  in 
Figure  26). 
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Haar  Wavelet 


Level 

Method  1 

Method  2 

Method  3 

Method  4 

3 

-90.34595401 

-70.62104470 

-80.60305640 

-124.69071864 

4 

-88.42929362 

-70.62113959 

-77.06934195 

-124.69071864 

5 

-87.46709193 

-70.62039117 

-74.23882826 

-124.69071864 

Daubechies  Order  3  Wavelet 


Level 

Method  1 

Method  2 

Method  3 

Method  4 

3 

-94.14439650 

-70.61933990 

-80.53273178 

-124.69071864 

4 

-92.47409703 

-70.61814619 

-76.96421916 

-124.69071864 

5 

-89.78405227 

-70.61533035 

-74.01710269 

-124.69071864 

Daubechies  Order  15  Wavelet 


Level 

Method  1 

Method  2 

Method  3 

Method  4 

3 

-96.53859183 

-70.61866118 

-80.44569092 

-124.69071864 

4 

-95.66683134 

-70.61739868 

-76.90756454 

-124.69071864 

5 

-91.83419067 

-70.60704310 

-73.94666010 

-124.69071864 

Coiflet  Order  1  Wavelet 


Level 

Method  1 

Method  2 

Method  3 

Method  4  ; 

3 

-92.87670574 

-70.61996912 

-80.48567609 

-124.69071864 

4 

-91.28425462 

-70.61972883 

-76.97295599 

-124.69071864 

5 

-88.52289832 

-70.60924047 

-74.05139371 

-124.69071864 

Meyer  Wavelet 


Level 

Method  1 

Method  2 

Method  3 

Method  4 

3 

-78.46926556 

-69.20423197 

-76.90477957 

-124.69071864 

4 

-75.28444411 

-68.17274988 

-73.30868711 

-124.69071864 

5 

-72.90028900 

-67.11307077 

-70.53191399 

-124.69071864 

Table  2.  Mean  Squared  Error  in  dB  Using  All  Three  Techniques  for  Various  Wavelet 
Basis  Functions  and  Multiple  Level  Decompositions  and  the  Fourier  method 

using  the  model  in  Figure  26. 


Similar  to  the  results  without  signal  extraction  (see  Table  1),  Method  4  produces 
the  most  similar  results  between  the  input  signal  and  desired  signal.  Additionally,  using 
Method  1  and  Daubechie’s  order  15  wavelet  with  three  levels  of  decomposition  provides 
the  best  results  when  compared  with  the  other  methods  and  wavelet  families.  In  this 
noise-free  environment,  the  signal  extraction  algorithm  perfonnance  is  slightly  improved 
from  the  non-extraction  case. 
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B. 


SIMULATION  MODEL  FOR  ANALYSIS  OF  SYSTEM  WITH  WHITE 
NOISE 


The  second  experiment,  based  on  the  models  depicted  in  Figures  27  and  28, 
captures  the  performance  of  both  the  wavelet  and  Fourier  methods  as  a  function  of  the 
signal-to-noise  ratios  (SNR);  the  SNR  is  the  ratio  of  signal  power  to  noise  power: 


SNR  = 


2>]2 

I>]2 


(5.1) 


where  x\n\  and  u{»]  are  the  signal  and  noise  samples,  respectively.  The  range  of  SNR 
evaluated  for  each  signal  was  from  -30  dB  to  30  dB,  but  for  each  plot,  the  results  are 
shown  only  for  the  range  of  interest. 


w[n\ 

Figure  27.  Basic  Model  of  System  Experiment  with  the  Addition  of  White  Noise. 

As  in  the  previous  section,  the  following  section  will  compare  results  with  and 
without  signal  extraction,  as  well  as  with  and  without  de-noising.  In  each  case,  the 
graphical  illustrations  capture  the  comparisons  between  Methods  1-3  and  Method  4.  The 
signal  and  channel  response  used  are  the  same  as  in  the  preceding  sections:  Waveform  1 
and  channel  response  “8694. cal”.  With  the  exception  of  adding  white  noise,  the  signal 
and  channel  response  preprocessing  steps  remain  the  same  as  described  above. 
Additionally,  the  mean  squared  error  in  each  plot  will  be  expressed  in  dB. 
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w[n\ 


Figure  28.  Basic  Model  of  System  Experiment  with  Additive  White  Noise  and 

Implementing  Signal  Extraction  Technique. 


1.  Equalization  Results  with  Additive  Noise  with  No  Signal  Extraction 

The  results  in  this  section  show  a  comparison  between  the  extracted  waveforms  as 
a  function  of  the  signal-to-noise  ratio  with  additive  white  noise.  It  uses  Waveform  1  and 
“8694. cal”  as  in  the  last  section,  as  well  as  the  same  algorithms  described  in  Chapter  IV. 
For  comparison  of  each  method,  Daubechie’s  Order  3  Wavelet  Level  3  was  used  for 
wavelet  decomposition  for  Method’s  1-3.  For  this  section,  the  signal  of  interest  is  not 
extracted  from  the  complete  recorded  signal. 

a.  Results  without  De-Noising 

In  this  case,  the  signal  has  not  been  extracted  and  de-noising  is  not  applied 
at  the  sub-band  level.  Among  all  four  methods,  there  is  little  difference  in  performance  as 
seen  in  Figure  29  (a). 
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Figure  29.  No  Signal  Extraction  and  De-noising.  (a)  Waveform  1  with  Additive 

Noise  Using  Methods  Using  Daubechie’s  Order  3  Wavelet  Level  3.  (b) 
Normalized  Results  Using  Daubechie’s  Order  3  Wavelet  Using  Level  3. 


Figure  29  (b)  shows  a  comparison  of  the  results  normalized  with  respect  to 
those  of  Method  4.  The  perfonnance  of  Methods  2  and  3  have  improved  over  that  of 
Method  4  for  SNR  less  than  3  dB  and  13  dB,  respectively.  Although  the  degradation  in 
performance  of  MSE  of  Methods  2-3  appears  to  exponentially  increase  above  the  0  dB 
point,  it  must  be  noted  that  the  scale  of  the  y-axis  is  on  the  order  of  thousandths  of  a  dB. 

b.  Results  using  De-noising 

When  the  same  experiment  was  conducted  for  a  signal  with  additive  noise 
and  no  signal  extraction,  but  with  de-noising  techniques,  the  results  are  significantly 
different.  As  can  be  seen  in  Figure  30  (a),  the  results  with  de-noising  have  notably 


49 


improved  for  SNR  less  than  17  dB.  However,  for  SNR  values  greater  than  17  dB,  the 
results  do  not  compare  closely.  These  results  are  not  a  significant  improvement  to  those 
above  where  de-noising  was  not  implemented  as  the  signal  to  noise  ration  of  the  signal 
increases. 


Figure  30.  No  Signal  Extraction,  (a)  Waveform  1  with  Additive  Noise  with  De- 

noising  Using  Daubechie’s  Order  3  Wavelet  Level  3.  (b)  Normalized  Results  with 
De-noising  Using  Daubechie’s  Order  3  Wavelet  Using  Level  3. 


The  time  domain  plots  for  this  case  are  displayed  in  Figure  3 1  at  different 
stages  of  the  processing.  It  illustrates  the  impact  of  additive  white  noise  and  visually 
shows  the  results  at  1 5  dB  SNR.  It  also  demonstrates  the  effect  with  and  without  de- 
noising  in  the  algorithm  between  the  input  signal,  s[n],  and  the  resultant  filtered 
signal,. ?[«]. 
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Figure  31.  Waveform  1  with  No  Signal  Extraction,  (a)  Original  Input  Signal,  (b) 

Signal  with  Noise  (SNR  of  15  dB).  (c)  Equalized  Signal  Using  Method  2  without 
De-noising  Using  Daubechie’s  Order  3  Wavelet,  (d)  Equalization  Using  Method  2 
with  De-noising  Using  Daubechie’s  Order  3  Wavelet. 


c.  Comparison  against  Other  Waveforms 

To  demonstrate  the  importance  of  selecting  a  wavelet  appropriate  for  the 
type  of  signal  being  processed,  the  following  comparison  shows  a  given  wavelet  against 
four  types  of  data.  Figure  32  shows  a  normalized  comparison  between  these  four  different 
waveforms  which  were  previously  shown  in  Figure  12.  The  experiment  was  conducted 
with  a  common  channel  response,  no  signal  extraction,  no  de-noising  and  using 
Daubechie’s  order  3,  and  3  levels  of  decomposition. 
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x  10  3 


Figure  32.  Comparison  between  Each  of  the  Waveforms  with  No  Signal  Extraction 

and  Normalized  MSE  Values  and  Using  Daubechie’s  Order  3,  and  3  Levels  of 

Decomposition. 


Figure  32  clearly  shows  the  impact  of  different  signals  with  a  given 
wavelet  family.  There  is  up  to  10  dB  of  degradation  in  performance  using  the  same 
wavelet  family,  but  applied  to  different  waveforms.  This  observation  indicates  that  better 
performance  can  be  achieved  by  selecting  wavelets  that  correlate  well  to  the  waveforms 
being  processed. 

2.  Equalization  Results  with  Additive  Noise  with  Signal  Extraction 

The  same  experiment  as  above  is  repeated  with  signal  extraction  as  depicted  in 
Figure  28,  using  Wavefonn  1,  “8694. cal”  as  a  fdter  response,  and  Daubechie’s  order  3, 
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and  3  levels  of  decomposition  for  Methods  1-3.  The  only  difference  is  the  use  of  the 
signal  extraction  method  discussed  in  Chapter  IV.  Effects  of  the  use  of  de -noising  are 
also  explored. 


a.  Results  without  De-Noising 

Figure  33(a)  shows  that,  using  the  configuration  noted  above  but  without 
de-noising,  there  is  limited  overall  improvement  as  the  SNR  is  increased.  A  closer  look  at 
Figure  33(b),  which  is  a  normalized  plot  with  respect  to  Method  4,  it  can  be  seen  that 
there  is  a  marked  improvement  of  performance  for  SNR  less  than  13  dB,  but  degradation 
occurs  at  high  SNR  values  (above  15  dB).  The  results  for  Methods  1  and  3  are  within  a 
hundredth  of  a  dB  for  this  case. 
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with  Additive  Noise  Using  Daubechie’s  Order  3  Wavelet,  Fevel  3.  (b) 
Normalized  Results  with  No  De-noising  Using  Daubechie’s  Order  3  Wavelet, 

Fevel  3. 
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b.  Results  with  De-Noising 

Figure  34  shows  results  with  de-noising  and  signal  extraction  using 
Waveform  1  and  “8694. cal”.  It  can  be  seen  again  that  there  is  improvement  in  the 
performance  for  SNR  less  than  approximately  15dB.  Compared  to  the  results  in  the 
previous  section  (see  Figure  33),  there  was  slight  improvement  to  the  performance  of  the 
system  at  higher  SNR  values  when  de-noising  techniques  were  implemented  as  shown  in 
Figure  34. 


Signal  Power  to  Noise  Power  (SNR)  in  dB 
(b) 


Figure  34.  MSE  as  a  function  of  SNR  with  Signal  Extraction  and  De-noising:  (a) 

Wavefonnl  with  Additive  Noise  Using  Daubechie’s  Order  3  Wavelet,  Level  3. 
(b)  Normalized  values  Using  Daubechie’s  Order  3  Wavelet,  Level  3. 
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A  last  comparison,  in  Figure  35,  shows  the  time  domain  representations  of 
the  waveforms  with  a  15  dB  SNR  using  extraction  techniques.  The  desired  signal  has 
been  extracted  from  the  data  containing  only  the  portion  of  the  data  where  the  signal  is 
present.  The  simulation  parameters  used  to  obtain  the  results  in  Figure  35  are  same  as 
those  used  to  obtain  the  results  in  Figure  3 1  with  the  exception  that  signal  extraction  is 
applied  here.  As  can  be  seen,  the  waveforms  have  been  processed  only  over  the  extracted 
portions. 


(a)  x10-s  (b) 


Time  in  sec  xIO-5  Time  in  sec  xIO"5 

(c)  (d) 

Figure  35.  Waveform  1  with  Signal  Extraction:  (a)  Original  Input  Signal,  (b)  Signal 

with  Noise  (SNR  of  15  dB).  (c)  Equalized  Signal  Using  Method  2  without  De¬ 
noting  Using  Daubechie’s  Order  3  Wavelet,  (d)  Equalization  Using  Method  2 
with  De-noising  Using  Daubechie’s  Order  3  Wavelet. 
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This  chapter  provided  results  and  plots  demonstrating  the  perfonnance  of 
the  algorithms  described  in  the  previous  chapters.  The  results  compared  a  single 
waveform  data  set  using  different  wavelet  types.  The  perfonnance  of  the  four  methods 
are  compared  as  a  function  of  the  signal-to-noise  ratio  (-30  to  30  dB)  using  wavelet  de¬ 
noting  techniques.  The  experiments  compared  the  difference  in  performance  with  and 
without  signal  extraction.  The  conclusions  of  the  work  reported  here  and  the  results  are 
discussed  in  Chapter  VI. 
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VI.  CONCLUSIONS 


This  thesis  applied  wavelet-based  signal  processing  techniques  to  an  existing 
practical  problem  of  aircraft  testing  for  vulnerability  to  EMP.  It  sought  to  investigate 
whether  the  performance  of  an  equalizer  could  be  improved  using  wavelet  decomposition 
versus  Fourier  transform,  whether  sub-band  de-noising,  as  well  as  preprocessing 
techniques  for  extracting  the  signal  from  arbitrary  data  would  prove  beneficial. 

A.  SIGNIFICANT  RESULTS 

First,  we  analyzed  whether  equalization  using  wavelet  transfonns  could 
outperform  equalization  using  Fourier  transforms.  The  investigation  included  five 
different  types  of  wavelets  and  three  different  levels  of  decomposition  each.  Results 
showed  that  Fourier  transforms  produced  a  closer  approximation  to  the  input  signal  after 
equalization.  This  can  be  attributed  to  the  fact  that  the  data  collected  from  the  channel 
was  already  in  the  frequency  domain,  and  the  scarcity  of  channel  response  data  and 
channel  dynamics  precluded  a  dramatic  increase  in  perfonnance.  Although  sub-band 
level  processing  is  expected  to  have  better  results,  given  the  known  channel  information, 
Fourier  transform  processing  was  less  computationally  intensive,  and  provided  a  better 
approximation  in  low  noise  environments. 

Second,  we  implemented  a  preprocessing  algorithm  to  remove  undesired  portions 
of  the  signal  from  a  given  signal  vector  using  time  average  energy  calculations.  Results 
showed  performance  improvement  and  a  closer  approximation  of  the  desired  signal  using 
this  technique.  Additionally,  from  a  test  and  evaluation  standpoint  in  order  to  support 
acceptance  and  certification,  it  will  keep  the  signal  quality  from  being  sensitive  to 
operator  intervention  and  error. 

Lastly,  we  investigated  the  equalization,  signal  extraction,  and  de-noising  for  a 
given  set  of  recorded  signal  measurements  with  the  addition  of  additive  noise.  Results 
showed  that  using  wavelets  had  better  performance  than  processing  in  the  frequency 
domain  for  signals  with  additive  noise;  however,  the  results  of  de-noising  are  preliminary 
and  further  in-depth  investigation  is  warranted  to  draw  definitive  conclusions  on  the 


57 


performance  improvement  with  de-noising.  For  each  algorithm  with  de-noising,  the 
wavelet  based  processing  had  better  results  for  signal-to-noise  ratios  of  less  than  15  dB. 

Results  suggest  that  the  current  method  employed  by  NAWCAD  for  channel 
equalization  using  Fourier  transforms  is  optimal  in  a  low  noise  (high  SNR)  environment 
based  on  the  channel  (equipment)  response  provided.  Although  potentially  beneficial 
with  other  systems,  wavelet  processing  did  not  provide  significant  improvement  to  the 
process  at  high  SNRs.  Additionally,  it  added  considerable  processing  overhead  for  the 
equalization  process.  Nevertheless,  significant  improvement  in  perfonnance  of  the 
proposed  wavelet  methods  with  de-noising  is  obtained  at  low  SNRs,  a  desirable  outcome 
especially  when  the  aircraft  testing  is  conducted  under  noisy  conditions.  Results  also 
suggest  that  removing  undesired  portions  of  a  signal  from  a  given  data  set  as  part  of 
preprocessing  provided  improved  system  performance.  By  automating  the  signal 
extraction  process  using  the  principle  of  time  average  energy,  one  can  obtain  consistency 
and  repeatability  of  the  results. 

B.  RECOMMENDATIONS  FOR  FURTHER  RESEARCH 

In  this  thesis,  using  basic  wavelet  de-noising  techniques  there  was  an 
improvement  in  signal  extraction  performance  in  high  noise  environments.  However,  the 
techniques  used  were  basic  and  no  in-depth  study  was  conducted.  A  more  detailed  look 
into  which  thresholding  methods  would  be  optimum  for  these  processes  will  be  useful. 
More  detailed  investigation  into  which  wavelet  families  work  best  for  these  classes  of 
signal  and  channel  structures  could  also  be  conducted  to  improve  system  perfonnance. 

The  ability  to  exploit  the  multi-rate  sampling  techniques  of  the  channel 
characteristics  is  another  viable  option.  Because  of  the  compression  capabilities  of 
wavelets,  although  not  discussed,  a  multi-rate  technique  may  be  developed  that  could 
capture  more  information  about  channel  using  fewer  points. 

Finally,  the  use  of  fractional  Fourier  transforms,  or  other  linear  or  non-linear 
transfonns,  may  also  offer  an  alternative  method  for  improving  the  performance  of  this 
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system.  Based  on  the  uniqueness  of  the  system  and  signal  characteristics,  the  solution 
should  not  be  limited  to  wavelet  or  Fourier  transform  processing  to  provide  the  most 
optimal  results. 
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APPENDIX 


A.  MAIN  PROGRAMS 

The  following  section  provides  the  MATLAB  code  used  for  the  results  generated 
in  Chapter  V.  It  provides  the  code  both  for  the  Tables  and  plots.  The  succeeding  section, 
Section  B,  provides  the  code  for  the  subroutines  called  in  the  following  program  sets. 

1.  Program  Producing  Tables  for  Chapter  V 


%  This  MATLAB  file  is  the  main  operating  program  for  simulations  of 
%  signal  processing  for  this  thesis. 

%  It  is  the  program  that  implements  the  extraction  technique  on  data 
%  and  presents  the  data  with  and  without  de-noising  techniques. 

%  The  results  of  this  program  produce  the  Tables  of  Chapter  V,  and  only 
%  results  for  the  original  signal  with  no  added  noise.  This  program 
%  primarily  compares  the  differen  t  wavelet  basis  functions  against 
%  one  another  as  well  as  compares  the  effectiveness  of  different  levels 
%  of  decomposition 
% 

%  Functions  Called: 

%  calprog2() 

%  buildiir() 

%  avepowerQ 

%  filterNreconstruct() 

%  polywavfill() 

%  wavedecom() 

%  noisepolywavfil  1  () 

%  noisewavedecom() 

% 


%  Clear  variables  and  close  all  figures 
clear;  close  all; 

%  Define  datal  as  global 
global  datal 


%  Loads  in  data 
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loadprogdatahome; 

%  Find  length  of  data  signal  to  be  used 
n = length  (data  1  (:,  2)); 

%o  Finds  the  freq  response  of  channel  response  data  @  freq  bins  of  input  data 
[caldata]  =calprog2(y4,deItafl  ,n); 

%  Finds  IIR  coefficients  and  HR  filter  frequency  response  of  order  40  from 

%  the  channel  response  data 

[g,  b  Hr,  a  Hr] = b  u  il  diir( ca  l da  ta(  I :  n/2 + 1 40,  ’n); 

%o  Established  full  spectrum  (both  halves  of  butterfly) 
gl  =[g '  fliplr(conj(g  (2:  length  (g)-l) '))] 

%o  Establishes  noise  contribution  from  the  signal 
sigpower =var( datal(:,2))  ; 

%  Establish  the  noise  power  to  be  added  to  the  signal 
sigma=0; 

%o  Creates  Data—  The  vector  to  be  filtered 
data=datal  ; 

%  Seeds  the  random  generator  to  keep  results  consistent  for  comparison 
randnC state' ,  I); 

%o  Adds  noise  to  signal 

data(:,2)=datal(:,2)  +  sigma*randn(length(datal(:,2)),l); 

%  Finds  the  extraction  parameters  using  the  function  avepower 
[dd,  start,  stop,  data ] =avepower( data); 

%  Extracts  the  signal  from  the  data  and  fills  in  zeros  to  match  the 
%  length  of  the  orignal  data 

datal(:,2)  =[datal  (start: stop, 2);  zeros (n-(stop-start+ 1) ,  1) ] ; 

%o  Plots  input  data  in  the  Time  and  Frequency  Domain 
figure('name',  'Original  Data-  Time  &  Freq  Domain'); 
subplot(2 ,1 ,1)  ;plot(datal  (:  ,1) ,  datal  (:  ,2),datal  (: ,  l),data(:  ,2)) 
subplot(2 , 1 ,2)  ;plot(fl  ,abs(fft(data  1  (:  ,2)))  fl  ,abs(fft(data(:  ,2)))) 

%  Creates  signal  x[n ]  that  needs  to  be  equalized 
%  Sends  input  throught  the  channel 
filterout=filter(biir,aiir,data(:  ,2))  ;disp('end  filter); 
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%  Equalizes  the  signal  x[n]  through  division  of  the  channel  frequency 
%  response.  Method  4  (their  method) 
ffout=fft(filterout) ; 
fyout=ffout./gl; 

%yout  is  shat[n]  for  method  4 
yout=reaI( iff't(fyout)); 

%  Plot  of  x[nj  in  Time  Domain 
figurefname',  'Filtered  Data  using  HR  Filter ) 
suhplot(2 , 1 , 1) ;plot(filterout) ;  axis  tight 

%  Plot  of  x[n]  in  Frequency  Domain 
subplot(2,l,2) 

plot((abs(ffout)));axis  tight;  title({ [Figure  1] .['filtered  output']}); 

%  Plot  of  shat[n]  in  Frequency  Domain 
figure  ('name',  'Equalized  Output  using  FFT) 

plot(fl,abs(fyout)fl,abs(fdatal(:,2)));  title({ ['Figure  1] ,[' equalized  output  usinf  [ft]}); 

%o  Initializes  parameters  for  wavelet  basis  functions  and  levels  of  decomposition 
wavnamelist={'haar'  'db3'  'dbl5'  'coifl'  'dmey'}; 
lev=[3  4  5] ; 

for  i=l:5 

wavname=wavnamelistf (}; 
disp('  );disp(' '); 

disp(sprintf('%-l  Os%',  wavname)) 

disp(sprintf('%-5s  %-15s  %-15s  %-15s  %-15s  %-15s  %-15s  %-l 5 s', 'level', 'Method 
1',  'Method  2',  'Method  3',  'Present  Tech’,  'Method  1  noise',  'Method  2  noise’,  'Method  3 
Noise')) 

for  j~  1:3 
level=lev(j); 

%  Find  Time  Domain  data  for  Approx' s  and  Details  for  current 
%  wavelet  basis  functions  and  number  of  levels 
[A,D]  =wavedecom(wavname, level Jilterout,  'n'j; 

%  Find  Time  Domain  data  for  Approx' s  and  Details  for  current 
%  wavelet  basis  functions  and  number  of  levels  using  de-noising 
%  techniques 

[ nA, nD]=noisewavedecom(wavname, level, filterout,  'n '); 

%o  Equalize  the  data  using  Methods  1  and  2,  for  the  with  and 
%  without  de-noising  methods. 

%  NEWDATA  is  the  equalized  output  using  inverse  HR  filter 
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%  (method  1) 

%  FFTNEWDATA  is  equalized  output  using  inverse  channel 
%  frequency  response  (method  2) 

[ n  e  wda  ta ,  ff'tn  e  wda  t  a ]  =filterNreconstruct(gl ,  aiir,  biir,A,D,  wavname); 
[noisenewdata.noisefftnewdata]  =filterNreconstruct(gl  ,aiir,biir,nA,nD,wavname) ; 

%  Equalization  of  the  data  Using  Method  3  for  the  with  and 
%  without  de-noising  methods. 

[polynewdatal ]  =polywavfill  (aiir,  biir,  wavname,  level, filter  out); 
[noisepolynewdatal]=noisepolywavfill  (aiir, biir, wavname, level, filter  out); 


%o  Find  the  Mean  Squared  Error  of  Techniques  1-3  (with  and  without  de-noising) 

%  in  dB,  using  original  noiseless  data  as  comparison 

si  =10*1  ogl  0(m  ean  ((datal(:,2) -new  data  ).A2)); 

s2=10*logl0(mean((datal(:,2)-Jftnewdata).A2)); 

s3=10*Iogl0(mean((datal(:,2)-polynewdatal).A2)); 

s4=10*logl0(mean((datal(:,2)-yout).A2)); 

s5=10*logl0(mean((datal(:,2)-noisenewdata).A2)); 

s6=10*logl0(mean((datal(:,2)-noisefftnewdata).A2 )); 

s7=10*logl0(mean((datal(:,2)-noisepoIynewdatal).A2)); 

disp(sprintf('%-5. 0f%-15.8f%-15.8f%-15.8f%-15.8f%-15.8f%-15.8f%-15.8f 
1 ,  level,  si  ,s2,s3  ,s4,s5  ,s6,s7)) 
end 
end 


2.  Program  that  Produces  Figures  for  Chapter  V 


%  This  MATLAB  file  is  the  main  operating  program  for  simulations  of 
%  signal  processing  for  this  thesis. 

%  It  is  the  program  that  implements  the  extraction  technique  on  data 
%  and  presents  the  data  with  and  without  de-noising  techniques. 

%  The  results  of  this  program  produce  the  Figures  of  Chapter  V,  and  only 
%  results  for  the  original  signal  with  added  noise  over  the  range  of  signal- 
%o  to-noise  ratio  (SNR)  from  -15  to  30  dB.  This  program  primarily  compares  the 
%  different  methods  against  one  another  for  a  given  wavelet  basis  function 
%  and  fixed  number  of  levels  for  wavelet  decomposition 
% 

% 

%  Functions  Called: 
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%  calprog2() 

%  buildiirQ 

%  avepowerQ 

%  filterNreconstruct() 

%  polywavfil  1  () 

%  wavedecom() 

%  noisepolywayfil  1  () 

%  noisewavedecom() 

% 


%  Clear  variables  and  close  all  figures 
clear;  close  all; 

%  Define  datal  as  global 
global  datal 

%  Loads  in  data 
loadprogdatahome; 

%  Find  length  of  data  signal  to  be  used 
n = length  (datal  (:,  2)); 

%  Finds  the  freq  response  of  channel  response  data  @  freq  bins  of  input  data 
[caIdata]=calprog2(y4,deItafl,n); 

%  Finds  IIR  coefficients  and  HR  filter  frequency  response  of  order  40  from 

%  the  channel  response  data 

[g,  biir,aiir] =buildiir(caldata(l  :n/2+l,  :),40,  'n '); 

%o  Established  full  spectrum  (both  halves  of  butterfly) 
gl  =[g ' fliplr(conj(g  (2:  length  (g)-l) '))]  '; 

%  Conduct  sweep  of  techniques  over  SNR  -15  to  30dB 
iter=30; 

snr=linspace(-15,30,iter); 

%  Initialize  Mean  Squared  error  matrices  to  save  results 
mmse  =zeros  (  7,  i  ter); 

%  Performs  the  channel  filtering  and  equalization  over  the  range  of  SNR  from  -15  to 
Ho  30  dB. 
for  v=l  liter 

%  Loads  in  original  data 
loadprogdatahome; 
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n=length(datal  (:,2)); 
close  all; 

%  Establishes  noise  contribution  from  the  signal 
sigpower=var(datal(:,2)); 

%  Establish  the  noise  power  to  be  added  to  the  signal 
sigma  =sqrt(sigpo  wer/ 1 0A(.l*snr(v))); 
disp(['The  amount  of  noise  is:  ’ ,num2str  (sigma)]); 

%  Creates  Data—  The  vector  to  be  filtered 
data=datal  ; 

%  Seeds  the  random  generator  to  keep  results  consistent  fot  comparison 
randnf state', 1); 

%  Adds  noise  to  signal 

data(:,2)  =datal(:,2)  +  sigma  *randn(length(datal  (:,2)),1); 

%  Finds  the  extraction  parameters  using  the  function  avepower 
[ dd,  start,  stop,  data] = avepower  (data); 

%  Extracts  the  signal  from  the  data  and  fills  in  zeros  to  match  the 
%  length  of  the  orignal  data 

datal  (: ,2)= [datal  (start: stop ,2) ;  zeros (n-(stop-start+ 1),  1 )] ; 

%  Plots  input  data  in  the  Time  and  Frequency  Domain 
figurefname',  'Original  Data-  Time  &  Freq  Domain '); 
subplot(2 ,1 ,1)  ;plot(datal  (:  ,1)  ,datal  (:  ,2)  ,datal  (:  ,1)  ,data(:  ,2)) 
subplot(2,l,2);plot(fl,abs(fft(datal(:,2))),fl,abs(fft(data(:,2)))) 

%  Creates  signal  x[n]  that  needs  to  be  equalized 
%  Sends  input  throught  the  channel 
filterout=filter(biir,aiir,data(:,2)); 

%o  Equalizes  the  signal  x[n ]  through  division  of  the  channel  frequency 
%  response.  Method  4  (their  method) 
ffout=fft(filterout); 
fyout=ffout./gl; 

%yout  is  shat[n]  for  method  4 
yout=real(ifft(fyout)) ; 

%o  Plot  of  xfn]  in  Time  Domain 
figurefname' ,  'Filtered  Data  using  HR  Filter) 
subplot(2 , 1 , 1) ;plot(filterout) ;  axis  tight 

%  Plot  of  xfn]  in  Frequency  Domain 
subplot(2,l,2) 

plot((abs(ffout)));axis  tight;  title({ [Figure  1] , [filtered  output']}); 
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%  Plot  of  shat[n]  in  Frequency  Domain 
figure('name',  'Equalized  Output  using  FFT) 

plot(fl,abs(fyout)fl,abs(fdatal(:,2)));  title({ [Figure  F] ,[' equalized  output  usinf  [ft]}); 

%o  Initializes  parameters  for  wavelet  basis  functions  and  levels  of  decomposition 
wavnamelist={'haar'  'db3'  'dbl5'  'coifF  'dmey'}; 
lev=[3  4  5] ; 

%o  Initialize  index  to  choose  wavelet  basis  function 
i=2; 

wavname=wavnamelistf i}; 
disp('  );disp(' '); 

disp(sprintf('%-l  Os%',  wavname)) 

disp(sprintf('%-5s  %-15s  %-15s  %-15s  %-15s  %-15s  %-15s  %- 15 s', 'level', 'Method 
F,  'Method  2’,  'Method  3\  'Present  Tech’,  'Method  1  noise',  'Method  2  noise',  'Method  3 
Noise')) 

%  Initialize  index  to  choose  number  of  levels  of  decomposition 

j=2; 

Ievel=lev(j); 

%  Find  Time  Domain  data  for  Approx' s  and  Details  for  current 
%  wavelet  basis  functions  and  number  of  levels 
[A,D]  =wavedecom(wavname,  level, filter  out,  'n '); 

%  Find  Time  Domain  data  for  Approx' s  and  Details  for  current 
%  wavelet  basis  functions  and  number  of  levels  using  de-noising 
%  techniques 

[ nA,  nDJ  =noisewavedecom(wavname,  level, filter  out,  ’n); 

%  Equalize  the  data  using  Methods  1  and  2,  for  the  with  and 
%  without  de-noising  methods. 

%  NEWDATA  is  the  equalized  output  using  inverse  HR  filter 
%  (method  1) 

%  FFTNEWDATA  is  equalized  output  using  inverse  channel 
%  frequency  response  (method  2) 

[newdatafftnewdata]  =filterNreconstruct(gl  ,aiir,biir, A, D,  wavname); 
[noisenewdata,noisefftnewdata]=filterNreconstruct(gl,aiir,biir,nA,nD, wavname); 

%  Equalization  of  the  data  Using  Method  3  for  the  with  and 
%  without  de-noising  methods. 

[polynewdatal]  =polywavfill  (aiir,biir, wavname, level, filterout) ; 

[noisepolynewdatal]  =noisepolywavfill  (aiir,biir, wavname, level, filterout) ; 

%o  Find  the  Mean  Squared  Error  of  Techniques  1-3  (with  and  without  de-noising) 

%  in  dB,  using  original  noiseless  data  as  comparison 
si  =10  *logl  0( mean(  (datal  (:,2)-newdata  ).A2)); 
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s2= 10*1  og  1 0(m  ean(( da  ta  1  (: ,  2)  -jftn  e  wda  ta ) .  '2 ) ) ; 
s3=l  0  *1 og  10( m  ean  ( ( da  ta  I  (: ,  2)  -po  lyn  e  wda  tal).  f2 ) ) ; 
s4=10  *logl  0(mean((datal  (:,2)-yout).  A2j); 
s  5=10*1  og  1 0(m  ean(( da  ta  I  (: ,  2)  -n  o  is  en  e  wda  ta).  f'2 ) ) ; 
s6=  10*1  og  1 0(m  ean(( da  ta  1  (: ,  2)  -no  is  efftn  e  wda  ta ) .  A2 ) ) ; 
s7=l  0  */ og  1 0(m  ean(( da  ta  I  (:,  2) -n  o  is  epo  lyn  e  wda  tal).  /'2) ) ; 

%o  Display  results  to  the  screen 

d  isp  (sprintf('%)-5. 0f%-15.8f%-15.8f%-15.8f%-15.8f%-15.8f%-15.8f%-15.8f 
'  level,  si,  s2,s3,s4,s5,s6,s7)) 

%  Save  current  iterations  data  to  matrix 
mmse(:,v)=[s  1  s2  s3  s4  s5  s6  s7]'; 

end 

%  Plot  of  Mean  Squared  Errors  of  Methods  over  SNR  from  -15  to  30  dB, 

%  N on-normalized  and  not  using  de-noising  techniques 
figure;  subplot(2, 1,1) 

plot(snr,mmse(l,:),snr,mmse(2,:),  'o-',snr,mmse(3,:),  'd-',snr,mmse(4,:),  's-');  axis  tight 
legend('Method  1 Method  2',  'Method  3',  'Present  Technique') 
xlaheK  1  ['Signal  Power  to  Noise  Power  (SNR)  in  dB] ,['(a)']j) 
ylabel('MSE  in  dB) 

%  Plot  of  Mean  Squared  Errors  of  Methods  over  SNR  from  -15  to  30  dB, 

%  Normalized  and  not  using  de-noising  techniques 
suhplot(2,l,2) 

plot(snr,mmse(l,:)-mmse(4,:),snr,mmse(2,:)-mmse(4,:),  'o- \snr,mmse(3,  :)-mmse(4, :), 'd- 
', sn r, m ms e( 4, :)-m ms e( 4,:),  's- ') ;  axis  tight 
legend('Method  1 ',  'Method  2',  'Method  3',  'Present  Technique) 
xlabel({ ['Signal  Power  to  Noise  Power  (SNR)  in  dB] ,['(b)']}) 
ylabel('MSE  in  dB) 

%  Plot  of  Mean  Squared  Errors  of  Methods  over  SNR  from  -15  to  30  dB, 

%  N on-normalized  and  using  de-noising  techniques 
figure;  subplot  (2, 1,1) 

plot(snr,mmse(5,:),snr,mmse(6,:),  'o-' ,snr,mmse(7 'd-',snr,mmse(4,:),  's-'J;  axis  tight 
legend('Method  1 ',  'Method  2\  'Method  3',  'Present  Technique) 
xlabel({ ['Signal  Power  to  Noise  Power  (SNR)  in  dB] ,['(a)']}) 
ylabel('MSE  in  dB) 

%  Plot  of  Mean  Squared  Errors  of  Methods  over  SNR  from  -15  to  30  dB, 

%  Normalized  and  using  de-noising  techniques 
subplot(2,l,2) 

plot(snr,  mmse( 5,  :)-mmse( 4,  :),snr,mmse( 6,  :)-mmse( 4, :),  'o-  ',snr,  mmse(7,  :)-mmse( 4, :), 'd- 
',  snr,  mmse( 4,  :)-mmse( 4, :), 's- ');  axis  tight 
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legend('Method  1  ’,  'Method  2',  'Method  3',  'Present  Technique') 
xlabel({ ['Signal  Power  to  Noise  Power  (SNR)  in  dB]  ,['(b)']j) 
ylabeI('MSE  in  dB) 


B. 


FUNCTIONS  USED  IN  ABOVE  MAIN  PROGRAMS 


1. 


AVEPOWER 


function  [dedrift, startlndex,  stoplndex,  dataExtracted]  =avepower(data); 
%%%0%%%0%%0%0%%0%%%%%0%0%0%%%%%0%%0%%0%%%0%%0%0%%0%%%%%0%%% 
%%%%%%%%%%%%%%%%%%%%% 

%  AVEPOWER 
% 

%  This  function  takes  a  data  sequence  "data"  and  returns 
%  signal  extracted  version  of  data  signal 

as  well  as  the  index  where  the  signal  starts, 
the  index  where  the  signal  stops,  and  the  mean 
of  the  values  not  contained  in  the  signal 


% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%%%%%%%%%%%%%%%%%%%%%%%% 


Inputs: 

data  (data  signal  to  be  extracted) 

Outputs: 

dedrift  (DC  Bias  of  non-signal  data)) 

Startlndex  (index  of  original  data  that  signal  begins) 
Stoplndex  (index  of  original  data  that  signal  ends) 
dataExtracted  (Extracted  signal  from  data) 


%  Establish  level  of performance  of  the  Threshold. 
threshperf=l; 

%  Determine  length  of  the  data 
n = length  (data (:,!)) ; 

%o  Establish  length  of  the  internal  to  average  over 
in  terra  l = ro  un  d( n/ /  0); 

%  Initialize  variable  for  energy  values 
energy=zeros(n,l); 

%o  Find  the  average  enrgy  values  over  length  of  data, 
for  i=l:n 

window=(i-(interval-l)); 

if  (i- (interval- 1))<1,  window=l;  end 
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energy  (i)=sum(data(window:i,  2).  A2)/ interval ; 
end 

%  Plot  the  Original  data 
figure 

subplot(4,l,l) 

plot(data(:  ,1)  ,data(:  ,2)) ;  axis  tight 

xlabel('(a)) 

ylabel(  Amplitude) 

%  Plot  the  Average  energy  calculation  of  the  data 
subplot(4,l,2) 

plot(data(:,l), energy, axis  tight 

xlabeK'(bj) 

y  lab  el (A  mpl  i  tude  ) 

%  Establish  thresholds  for  the  data,  threshl  represents  the  threshold  at 
%  the  beginning  of  the  data,  thresh2  represents  the  data  at  the  end  of 
%  the  data  string 
thresh  1  =  threshperf*l .  5; 
thresh2=threshperf*2 ; 

%o  Find  the  Beginn  ing  of  the  signal;  where  the  energy  crosses  the  threshold 

startlndex=interval; 

count=l; 

while  (threshl  ^ energy (startlndex)  >  ener gy (count)) &(count<n) 
count=count+ 1 ; 
end 

%  Add  portion  of  data  to  prevent  unnecessary  loss  of  the  signal 
startIndex=(count-round(count/5)); 

%>  Find  the  end  the  signal  in  the  data;  where  the  end  of  the  signal  crosses 

%  the  threshold. 

stoplndex=n; 

count=n; 

while  (thresh2*energy(n)  >  energy(count))&(count>l) 
count=count-l  ; 
end 

%  Add  portion  of  data  to  prevent  unnecessary  loss  of  the  signal 
stoplndex=(count+round(interval/2)); 

%  Check  values  to  make  sure  they  are  valid 
if  stoplndex>=n,  stopIndex=n;  end 
if  startlndex<l,  startlndex=l;  end 
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if  stopIndex<startIndex,  stopIndex=n;  startlndex=l;  end 

%  Eliminate  portions  in  the  average  time  energy  profile  that  do  not  con  tain  the  signal 
energy= [energy (startlndex: stoplndex) ;  zeros (n-(stopIndex-startIndex+ 1) ,1)] ; 

%  Eliminate  portions  of  the  data  that  does  not  contain  portions  of  the  signal 
dataExtracted(:  ,2)=  [data(startlndex: stoplndex, 2);  zeros(n-(stopIndex-startIndex+ 1),  1)J ; 

%  Calculate  the  average  noise  energy  of  the  data  not  considered  the  signal 
dcdrift=mean( [data (interval: startlndex, 2)  ;data(stopIndex:n,2)]); 

%  Plot  the  extracted  time  average  energy  of  the  signal 
subplot(4,l,3) 

plot (data(:,l), energy);  axis  tight 

xlabel('(c)) 

ylabel  (Amplitude) 

%  Plot  the  extracted  signal 
subplot(4,l,4) 

plot(data(:,l),dataExtracted(:,2));  axis  tight 
xlabel({['Time  in  sec] ,['(d)]}) 
ylabel  ('A  mpl  i  tude  ) 


2.  CALPROG2 


function  [calnew]  =calprog2 (cal,df,nx) ; 


%  CALPR0G2 

% 

%  This  function  takes  given  channel  response, 

%  This  code  was  modified  from  foldprob.m  [ ref] 

%  Inputs: 

%  cal  (data  signal  to  be  extracted) 

%  df  (delta  frequency  of  data  signal) 

%  nx  (number  of  data  points  in  the  data  signal) 

%  Outputs: 

%  calnew  (Extracted  signal  from  data) 

%  Function  calls: 

%  interpjkQ 

%  extendQ 

%  unwrap() 

% 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%%%%%%%%%%%%%%%%%%%%%%%% 
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%  INTERPJK  Table  look-up. 

%  Y  =  interpjk(datasetl,dataset2)  returns  a  dataset  of  linearly 

%  interpolated  rows,  mapping  dataset  1  to  the  first  column  of 
%  dataset2. 

% 


%  NOTE:  dataset2’s  1st  column  and  dataset  1  must  be  monotonically 
%  increasing. 


%  EXTEND 

%  dataset  1  =extend(datasetl ,dataset2) 

% 


%  Use  extend  to  extend  the  data  to  allow  for  interpolation  of 
%  datasets  whose  points  do  not  match  up. 


%  UNWRAP  Unwrap  phase  angle. 

%  UNWRAP  (P)  unwraps  radian  phases  P  by  changing  absolute 
%  jumps  greater  than  pi  to  their  2  *pi  complement.  It  unwraps 
%  along  the  first  non-singleton  dimension  of  P.  P  can  be  a 
%  scalar,  vector,  matrix,  or  N-D  array. 

% 

%  UNWRAP (P,  TOL)  uses  a  jump  tolerance  of  TOL  rather 
%o  than  the  default  TOL  =  pi. 

%0 


%  UNWRAP (P,  [], DIM)  unwraps  along  dimension  DIM  using  the 
%  default  tolerance.  UNWRAP(P,TOL,DIM)  uses  a  jump  tolerance 
%  of  TOL. 


%  Find  the  number  of  data  points  in  the  channel  frequency  response 
n=length(cal); 

%o  Set  up  frequency  bins  for  channel  response  interpolation 
f=linspace( 0,  df*nx/2,  nx/2 +1); 

%  Interpolates  data  and  sets  data  to  the  scale  of  the  vector  f 
%  Code  taken  from  foldprob.n  [ref] 
newdata=interpjk(cal,J ); 
cal_id=extend(cal,f; 

t= [cal _id(: ,  1)  ,abs(cal _id(:  ,2)) , angle  (cal _id( :  ,2))J  ; 
t(:,  3 )  =unwrap( t(:,  3)); 
ct=interpjk( t,f); 
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clear  t; 

t=ct(:,2).  *(cos(ct(: ,3))+sin(ct(: ,3))  *i); 
clear  caljd; 

%o  saves  modified  interpolated  data  as  newdata 
newdata=t'; 

%  Takes  the  one-sidedspectrum  and  adds  its  complex  conjugate  other  side 
if  length(newdata)/2  ~=  round(length(newdata)/2) 
fulldata= [newdata  fliplr(conj  (newdata(2 :  length  (newdata)- 1)))] ; 
else 

fulldata=  [newdata  fliplr(conj(newdata(2 :  length(newdata))))]  ; 
end 

%  Saves  final  channel  response  as  calnew 

calnew=[linspace(0,(length(fulldata)-l)  *df,length(fulldata)) '  fulldata] ; 


%  Plot  of  channel  frequency  response  before  and  after  interpolation 
figure  ('name',  'Channel  Freq  Domain-  Before/  After  interpolation') 
subplot(2,l,l) 

plot(cal(:,l),abs(cal(:,2)), '.') 

axis  tight 

subplot(2,l,2) 

plot (f,abs (newdata), '. ) 

axis  tight 


%  Plot  of full  spectrum  of  channel  response  and  its  time  domain  impulse 
%  response. 

figure  ('name',  'Channel  Freq  and  Time  Domain) 
subplot(2,l,l) 

plot(calnew(:  ,1)  ,abs(calnew(:  ,2))) 
subplot(2,l,2) 

[f=ijft  (fulldata); 

plot(abs(ff)); 

axis([0  4000  0  .0001]) 


3.  WAVEDECOM 


function  [A, D]=wavedecom(wavlet,  level,  data, pr); 


% 

%  This  function  finds  the  time  domain  sub-band  signals,  Approximation  and 
%  Detail  for  a  given  data  vector  "data  "  with  "level"  number  of  levels 
%  decomposition  and  using  the  "wavelet"  basis  function.  All  Approximation 
%  and  Detail  sub-band  signals  are  found  for  all  levels  and  are  arranged 
%  in  "level"  number  of  vectors  in  the  A  and  D  Matrices,  respectively 
% 

%  Inputs: 

%  wavlet  (Type  of  Wavelet  basis  function  to  be  used) 

%  level  (Number  of  levels  of  decomposition) 

%  data  (Data  vector  to  be  decomposed) 

%  pr  (string  variable,  'y',  prints  all  associated  plots) 

%  Outputs: 

%  A  (sub-band  Time  Domain  Approximation  signals  in  matrix  form  of  all  levels) 

%  D  (sub-band  Time  Domain  Detail  signals  in  matrix  form  of  all  levels) 


% 


%  Determine  the  length  of  the  data 
n=length(data); 

%  Decompose  the  data  using  above  described  parameters 
[C,L]  =  wavedec(data, level, wavlet) ; 

%  Find  the  Approximation  coefficients  for  last  level  of  decomposition 
Atmp=  appcoef(C,L, wavlet, level); 

%o  Zero  pad  the  remaining  length  of  Approximation  coefficients  up  to  length  n 
n  1  =  1 :  n/(2 Alevel)  ; 

A  =[Atmp '  zeros  (1,  n-L(l))] ; 

%o  Find  the  Detail  coefficients  for  all  levels  of  decomposition 
Dc=detcoef(C,L,  wavlet,  level) 

%o  If  requested  Display  the  following  plots 
if  pr  ==  y 

%  Plot  of  the  last  level  of  Approximation  coefficients  and  all  levels 
%  of  Detail  Coefficients 
figurefname',  'Wavelet  Coefficients') 
subplot(2,2,l) 
plot(Atmp);axis  tight 
titIe(['A',num2str(level), '  Coefficients'] ) 

subplot(2,2,2) 
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plot(Dc{level});axis  tight 

title(['D \num2str (level), '  Coefficients'] ) 

subp!ot(2,2,3) 

pIot(Dc{leveI-lj);axis  tight 

title(['D[ ,num2str  (level- 1), '  Coefficients'] ) 

subplot(2,2,4) 

p\ot(Dc{level-2});axis  tight 

title( ['D' ,num2str (level-2) , '  Coefficients'] ) 

for  i=l:level-3 
if  mod(i-l,4)==0 

figure('name',  'Wavelet  Coefficients') 
end 

subplot(2,2,mod(i-l,4)+l);  plot(Dc{IeveI-(i+2)});  axis  tight 
titIe(['D',num2str(leveI-(i+2)), '  Coefficients'] ) 
end 

%  Plot  of  Original  Signal 
figure('name',  'Original  data) 
plot(data);axis  tight 
end 

%  Initialize  matrices 
A  =zeros(level,n); 

D  =zeros  (level,  n); 

%  For  each  level  of  decomposition  save  the  selective  reconstruction  of  the 
%  Approximation  and  Detail  Time  Domain  signals 
for  i=l:  level 

A(i, :)  =wrcoef('a ',  C,L,  wavlet,  i) '; 

D(i,:)  =  wrcoef('d',C,L, wavlet, i)'; 
end 

%  If  requested,  display,  the  following  plots  of  Time  domain  sub-band 
%  signals  for  all  Approximation  and  Detail  levels  of  decomposition 
if pr=='y' 
for  i=l  .level 

if  mod(i-l,4)==0 

figure('name',  'Wavelet  Time  Domain  Approximations  and  Details) 
end 

subplot(4,2,mod(2*(i-l),8)+l); plot(A(i,:));  axis  tight 
title(['A',num2str(i),'  Reconstruction  to  Time  Domain']) 

subplot(4,2,mod(2*(i-l),8)+2); pIot(D(i,:));  axis  tight 
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title(['D',num2str(i), '  Reconstruction  to  Time  Domain']) 
end 
end 


4.  BUILDIIR 


function  [g,b,a]  =buildiir(cal, order, pr) 


%  BUILDIIR 

% 

%  This  function  finds  HR  filter  coefficients  for  the  best  fit  filter 
%  given  frequency  response  using  a  mean  squared  error  criterion  . 

%  Yule-Walker  function  from  MATLAB  used  to  determine  coefficients  from 
%  Frequency  Response. 

% 

%  Inputs: 

%  cal  (FFT  of  channel  response) 

%  order  (Order  of  desired  HR  filter,  zero  implies  MMSE  best  fit) 

%  pr  (string  variable,  'y',  prints  all  associated  plots) 

%  Outputs: 

%  g  (Frequency  Response  of  HR  Filter) 

%  a  (HR  Filter  denominator  coefficients) 

%  b  (HR  Filter  numerator  coefficients) 

% 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%%%%%%%%%%%%%%%%%%%%%%%% 


%  Fimd  the  length  of  channel  frequency  response  vector 
n=length(cal(:,l)); 

%o  Normalize  the  Channel  Frequencies  from  0  to  1 
f=cal(:,l)/cal(n,l); 

%o  Establish  vestor  of  Channel  Magnitudes 
m=abs(cal(:,2)); 

%  Find  the  optimum  filter  if  the  input  order  is  zero 
if  order==0 

%  Initialize  the  comparison  of  the  error  value 
besterror=10Al  6; 

%  Find  the  best  filter  from  order  2  to  80,  with  even  order  values. 
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for  i=l:40 

%  Find  filter  coefficients  of  order  2  *i  using  yulewalk 
[b  I  ,a  I ]  =  yulewalk(2  *ifm); 

%o  Find  frequency  response  of  resultant  HR  filter  using  freqz 
[h,wj  =  freqz ( bl,a  I,  length(f) ) ; 

%o  Compare  the  difference  of  HR  filter  and  Channel  response  using 

%  MSE  comparison. 

error =m  can  ((a  bs  ( h)  -a  bs  ( m ) ) .  r'2 ) ; 

%o  If  the  response  was  better  than  order  2*(i-l)  then  save  the  order 
%  Also,  verify  that  roots  and  poles  of  HR  filter  are  within  unit 
%  circle  for  stability. 

if  (error  <  b  es  terror)  &(  (a  bs  ( ro  ots(a  I )))<I)&((a  bs  (roots(b  l )))<!) 
order=i*2; 
besterror=error; 
end 
end 
end 

%  Calculate  Filter  coefficients  fro  best  order  HR  Filter 
[bl  ,al]  =  yulewalk(orderfm); 

%o  Calculate  Frequency  Response  for  best  filter 
[h,wj  =freqz(bl,al,length(f)); 

%o  If  requested,  plot  the  following  results 
ifpr  ==  y 

%  Plot  the  Original  channel  response  and  the  HR  filter  frequency 
%  response 

figure('name ',  HR  Freq  Response') 
plot(f,abs(m),w/pi,abs(h), '—) 
legend(’ Original' ,  'yulewalk  Designed) 

title({ ['Figure  1 —  HR] ,  ['Comparison  of  Frequency  Response  Magnitudes'] }) 

%  Plot  the  zeros  and  poles  of  HR  Filter 
figure  ('name',  'IIR  Poles  and  Zeroes) 
zplane(bfal) 

title(’Poles  and  Zeros  of  IIR  Filter) 
end 

%  Returned  values 

b=bl; 

a=al; 

g=h; 
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5. 


FILTERNRECONSTRUCT 


function  [newdata,fftnewdata]=filterNreconstruct(fftfil,b,a,A,D,wavlet) 


%  FILTERNRECONSTRUCT 


% 

%  This  function  performs  Methods  1  and  2  on  the  data  for  Equalization 
% 

%  Inputs: 

%  fftfil  (FFT  channel  response) 

%  a  (HR  Filter  coefficients) 

%  b  (HR  Filter  coefficients) 

%  A  (Matrix  of  Approximation  Time  Domain  sub-band  Signals) 

%  D  (Matrix  of  Detail  Time  Domain  sub-band  Signals) 

%  wavelet  (wavelet  basis  function  to  be  used) 

%  Outputs: 

%  newdata  (Equalized  data  using  Method  1) 

%  fftnewdata  (Equalized  data  using  Method  2) 

% 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%%%%%%%%%%%%%%%%%%%%%%%% 


%  Determines  the  size  of  Matrix  A  to  determines  level  of  decomposition 
%  and  length  of  data  string 
[level,  n]  =size(A); 

%o  Filters  each  sub-band  signal 
for  i=l:  level 

%  Filters  each  sub-band  using  inverse  HR  filter 
filA(i,:)  =filter(b,a,A(i,:)) ; 
filD(i, :)  =  filter  (b,  a,D( % :)); 

%o  Filters  each  sub-band  using  the  inverse  of  the  frequency  response 
fftfilAfi :)  =real(ifft(fft(A(i, :). /fftfil '))); 
fftfilD(i, :)  =real(ifft(fft(D( i, :). /fftfil '))); 
end 


%  Decomposes  Each  sub-band  signal  into  its  respective  coefficients  for  HR 
%  method  (Method  1) 

[C,L]  =  wavedec(filA(level, :), level, wavlet) ; 

Atmp=  appcoef(C,L, wavlet, level); 

Cnew=Atmp; 
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%  Decomposes  Each  sub-band  signal  into  its  respective  coefficients  for 
%  frequency  response  method  (Method  2) 

[fC,fL]  =  wavedec(fftfilA(level, :),  level,  wavlet) ; 
fAtmp=  appcoef(fC,fL,  wavlet,  level); 
fCnew=fAtmp; 

%  Puts  together  coefficients  of  each  sub-band  to  reconstruct  signal 
for  i=l:  level 

%  Method  1  signal 

[C,L]  =  wavedec(filD(level-i+ 1 ,:),  level,  wavlet) ; 
tmp=  detcoef(C,L, wavlet, level); 

Cnew=[ Cnew  tmp{level-i+l}J; 

%  Method  2  signal 

[fC,fL]  =  wavedec(jftfilD(level-i+l,:),  level,  wavlet); 
ftmp  =  detcoef(fC,fL, wavlet, level); 
fCnew= [fCnew  ftmp{level-i +1} ]; 
end 

%  Reconstructs  the  complete  signal  from  coefficients 
newdata  =  waver  ec(  Cnew,L,  wavlet); 
fftnewdata  =  waverecf Cnew, fL, wavlet) '; 


6.  POL  YWAVEFIL 1 


function  [polynewdata]  =polywavfil  1  (a, b, wavlet, level, filterout) ; 


%  POLYWAVFIL1 

% 

%  This  function  performs  Method  3  on  the  data  for  Equalization 
% 


% 

% 

% 

% 

% 

% 

% 

% 

% 


Inputs: 

a  (IIR  Filter  coefficients) 
b  (IIR  Filter  coefficients) 
filterout  (signal  to  be  equalized) 
level  (number  of  levels  of  wavelet  decomposition) 
wavelet  (wavelet  basis  function  to  be  used) 
Outputs: 

polynewdata  (Equalized  data  using  Method  3) 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%%%%%%%%%%%%%%%%%%%%%%%% 


79 


%  Decomposes  the  data  signal  into  coefficients 
[pC,pL]  =  wavedecfiilterout,  level,  wavlet); 

%  Initializes  index  and  new  coefficient  vectors 
index=l; 

PolyC-pC;  PoIyL=pL; 

%  Filters  the  coefficients  of  the  data  signal  by  the  inverse  HR  filter 
for  i=l :  level +1 

tmpL=pL(i); 

tmp=(pC(index:index+tmpL-l)) ; 

PolyC(index:index+tmpL-l)=filter(a,b,  tmp); 
index=index+tmpL; 
end 

%  Reconstructs  the  resultant  data  using  new  filtered  coefficients 
polynewdata=  waverec(PolyC,PolyL, wavlet) 


C.  NOTES  ON  DE-NOISING  AND  EXTRACTION 
1.  De-Noising 


%Q%Q%0%Q%Q%Q%Q%0%0%Q%Q%Q%Q%0%Q%Q%Q%Q%0%0%Q%Q%0%Q%Q%Q%0%Q%Q%Q%0%0%Q%Q%Q%0%0%0%Q%Q%Q%0%Q 


%%%%%%%%%%%%%%%%%%%%% 

%  NOISEWA  VEDECOM  and  NOISEPOLYWA  VFIL1 


% 

%  These  two  functions  perform  the  same  as  wavedecom  and  polywavfil  1 , 
%  respectively  with  the  exception  of  one  line  of  code.  In  these 
%  functions,  when  the  data  is  decomposed  into  wavelet  coefficients, 

%  appropriate  de-noising  techniques  are  applied  as  well.  All  input 
%  and  output  parameters  remain  the  same.  The  modifications  of  each  of 
%  these  functions  are  show  below 
% 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%%%%%%%%%%%%%%%%%%%%%%%% 


% 

%  Lines  of  code  in  wavedecom 

%  Decompose  the  data  using  above  described  parameters 
[C,L]  =  wavedec  (data,  level,  wavlet); 

%  are  replaced  with  the  following  lines  of  code: 

%  Decompose  the  data  using  above  described  parameters  as  well  as 
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%  de-noising  techniques  with  parameters  'huersure'  and  's' 
[x,C,L] =wden(data,  'heursure',  's',  'one', level, wavlet); 


% 

%  Lines  of  code  in  polywavfill 

%  Decomposes  the  data  signal  into  coefficients 
[pC,pL]  =  wavedec(filterout, level, wavlet); 

%  are  replaced  with  the  following  lines  of  code  in  noisepolywavfill 
%  Decompose  the  data  using  above  described  parameters  as  well  as 
%  de-noising  techniques  with  parameters  'huersure'  and  's' 
[x,pC,pL]  =wden(filterout,  'heursure',  's',  'one' , level, wavlet) ; 


2.  Extraction 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%%%%%%%%%%%%%%%%%%%%% 


%  Extraction 

% 

%  The  operation  of  extracting  the  signal  from  noise  was  accomplished 
%  through  coding  with  the  addition  of  the  function  avepower  to  the 
%  program  code.  If  these  corresponding  lines  of  code  were  removed  and 
%  replaced  with  alternate  lines,  the  program  would  operate  without  the 
%  abilitv  to  extract  the  signal. 

% 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%%%%%%%%%%%%%%%%%%%%%%%% 


%  Lines  of  code  responsible  for  extraction  of  signal 

%  Finds  the  extraction  parameters  using  the  function  avepower 
[ dd,  start,  stop,  data] = avepower  (data); 

%o  Extracts  the  signal  from  the  data  and  fills  in  zeros  to  match  the 
%  length  of  the  original  data 

data  1  (:,  2)= [data  1  (start: stop,  2);  zeros (n-(stop-start+ 1),  1)] ; 

%  Lines  of  code  required  in  place  of  above  code 
datal=data; 
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